/*********************************************************************
*                                                                    *
* derivs_3d.c -- take derivatives of fields using second or fourth   *
*                order method                                        *
*                                                                    *
* Parameters:                                                        *
*        grid    -- pointer to grid values structure                 *
*       field    -- pointer to field to derivate                     *
*      dfield    -- pointer to field where derivative is to be stored*
*  coordinate    -- pointer to coordinate with respect to which      *
*                   the derivative is taken                          *
* Returns: nothing                                                   *
* Status: tested OK                                                  *
*                                                                    *
*********************************************************************/




/* #include "gen_3d.h" */
#include <math.h>
#include "first_macro_2d.h"
#include "structs_2d.h"    /* place where the used structures are defined */

#ifdef DERIVS_2D_H
#else
#define DERIVS_2D_H

/* SECOND ORDER ONES */

void derivS_2d(struct grid_2d *grid,
	     struct field *field,
	     struct field *dfield,
	     int *coordinate)
{/* derivS_2d */

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;

    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;

#ifdef EXCISION
    int exc_n1 = (*grid).exc_grid_1;
    int exc_n2 = (*grid).exc_grid_2;
#endif

    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;


    FLOAT dos = 2.0;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi)*0.5;
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi)*0.5;

#ifdef EXCISION
    int exc_n1i = (n1f-n1i-1)/2-exc_n1+1; /* take out a -1 for visual compatibility with the rest of code */
    int exc_n2i = (n2f-n2i-1)/2-exc_n2+1;
    int exc_n1f = (n1f-n1i-1)/2+exc_n1; 
    int exc_n2f = (n2f-n2i-1)/2+exc_n2;
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{   /*  case 1: */

	    {register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = dos*(-(*field).u_comp[IND_2(n1i,k)]
					           +(*field).u_comp[IND_2(n1i+1,k)]
						   )*h_1;


		/* Intermediate Values 1---(n-2) */

		{register int j;
		for (j = n1i+1; j < n1f-1; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (- (*field).u_comp[IND_2(j-1,k)]
					         + (*field).u_comp[IND_2(j+1,k)]
						 )*h_1;
		}
		}

		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] = dos*(+ (*field).u_comp[IND_2(n1f-1,k)]
					             - (*field).u_comp[IND_2(n1f-2,k)]
						     )*h_1;


	      }
	    }
	    

#ifdef EXCISION

	    /* interior and fases */
	  {register int n1;
	  register int n2;
	  for (n2 = exc_n2i; n2 < exc_n2f; n2++){
	      for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		(*dfield).u_comp[IND_2(n1,n2)] = 0.0;
	    }
	      (*dfield).u_comp[IND_2(exc_n1i-1,n2)] = dos*((*field).u_comp[IND_2(exc_n1i-1,n2)] 
		                                 - (*field).u_comp[IND_2(exc_n1i-2,n2)])*h_1;
	      (*dfield).u_comp[IND_2(exc_n1f,n2)] = dos*((*field).u_comp[IND_2(exc_n1f+1,n2)] 
		                                   - (*field).u_comp[IND_2(exc_n1f,n2)])*h_1;
	    }
	  }
	  

	    /* edges different from the one taking derivatives */

		(*dfield).u_comp[IND_2(exc_n1i-1,exc_n2i-1)] = - dos*(
                                 - (*field).u_comp[IND_2(exc_n1i+0,exc_n2i-1)] 
				 - (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-1)]        
			     + dos*(*field).u_comp[IND_2(exc_n1i-2,exc_n2i-1)]
				 )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1i-1,exc_n2f)] =  - dos*(
                                 - (*field).u_comp[IND_2(exc_n1i+0,exc_n2f)] 
				 - (*field).u_comp[IND_2(exc_n1i-1,exc_n2f)]        
			     + dos*(*field).u_comp[IND_2(exc_n1i-2,exc_n2f)]
				 )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2i-1)] =    dos*(
                               dos*(*field).u_comp[IND_2(exc_n1f+1,exc_n2i-1)] 
		                 - (*field).u_comp[IND_2(exc_n1f+0,exc_n2i-1)]        
				 - (*field).u_comp[IND_2(exc_n1f-1,exc_n2i-1)]
			       )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2f)] =    dos*(
                               dos*(*field).u_comp[IND_2(exc_n1f+1,exc_n2f)] 
			         - (*field).u_comp[IND_2(exc_n1f,exc_n2f)]        
			         - (*field).u_comp[IND_2(exc_n1f-1,exc_n2f)]
			       )*h_1/3.0;
	    


#endif	 



#ifdef DEBUG
	    printf("at case 1");
#endif
	    }

	break;


	/* Take derivative of second coordinate */

	case 2:{/* 	case 2: */


	    {register int k;
	    for (k = n1i; k < n1f; k++){
		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = dos*(- (*field).u_comp[IND_2(k,n2i)]
					           + (*field).u_comp[IND_2(k,n2i+1)]
						   )*h_2;



		/* Intermediate Values 1---(n-2) */

		{register int j;
		for (j = n2i+1; j < n2f-1; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (- (*field).u_comp[IND_2(k,j-1)]
					      	 + (*field).u_comp[IND_2(k,j+1)]
						 )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = dos*(+ (*field).u_comp[IND_2(k,n2f-1)]
					       	     - (*field).u_comp[IND_2(k,n2f-2)]
						     )*h_2;
	      }
	    }
	    

#ifdef EXCISION

	  {register int n1;
	  register int n2;
	    for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		for (n2 = exc_n2i; n2 < exc_n2f; n2++){
		  (*dfield).u_comp[IND_2(n1,n2)] = 0.0;
		}
		(*dfield).u_comp[IND_2(n1,exc_n2i-1)] = dos*((*field).u_comp[IND_2(n1,exc_n2i-1)] 
		                                   - (*field).u_comp[IND_2(n1,exc_n2i-2)])*h_2;
		(*dfield).u_comp[IND_2(n1,exc_n2f)] = dos*((*field).u_comp[IND_2(n1,exc_n2f+1)] 
						     - (*field).u_comp[IND_2(n1,exc_n2f)])*h_2;
	      }
	    }
	  

	    /* edges different from the one taking derivatives */



		(*dfield).u_comp[IND_2(exc_n1i-1,exc_n2i-1)] =  - dos*(
                                  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-0)] 
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-1)]        
				  + dos*(*field).u_comp[IND_2(exc_n1i-1,exc_n2i-2)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1i-1,exc_n2f)] =    dos*(
                                  + dos*(*field).u_comp[IND_2(exc_n1i-1,exc_n2f+1)] 
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2f+0)]        
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2f-1)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2i-1)] =  - dos*(
                                  -     (*field).u_comp[IND_2(exc_n1f,exc_n2i-0)] 
				  -     (*field).u_comp[IND_2(exc_n1f,exc_n2i-1)]        
				  + dos*(*field).u_comp[IND_2(exc_n1f,exc_n2i-2)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2f)] =    dos*(
                                    dos*(*field).u_comp[IND_2(exc_n1f,exc_n2f+1)] 
				    -  (*field).u_comp[IND_2(exc_n1f,exc_n2f)]        
				    -  (*field).u_comp[IND_2(exc_n1f,exc_n2f-1)]
				    )*h_2/3.0;
	    

  
#endif	 


#ifdef DEBUG
	    printf("at case 2");
#endif
	}
	
	break;


	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}

	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}

}




/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */



void derivS_Per_x_2d(struct grid_2d *grid,
	     struct field *field,
	     struct field *dfield,
	     int *coordinate)
{/* derivS_2d */

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;

    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;

#ifdef EXCISION
    int exc_n1 = (*grid).exc_grid_1;
    int exc_n2 = (*grid).exc_grid_2;
#endif

    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;


    FLOAT dos = 2.0;

    FLOAT h_1 = (FLOAT)(n1f-n1i)/(xf-xi)*0.5; // Periodic
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi)*0.5;

#ifdef EXCISION
    int exc_n1i = (n1f-n1i-1)/2-exc_n1+1; /* take out a -1 for visual compatibility with the rest of code */
    int exc_n2i = (n2f-n2i-1)/2-exc_n2+1;
    int exc_n1f = (n1f-n1i-1)/2+exc_n1; 
    int exc_n2f = (n2f-n2i-1)/2+exc_n2;
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{   /*  case 1: */

	    {register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = (-(*field).u_comp[IND_2(n1f-1,k)]
					           +(*field).u_comp[IND_2(n1i+1,k)]
						   )*h_1;


		/* Intermediate Values 1---(n-2) */

		{register int j;
		for (j = n1i+1; j < n1f-1; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (- (*field).u_comp[IND_2(j-1,k)]
					         + (*field).u_comp[IND_2(j+1,k)]
						 )*h_1;
		}
		}

		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] =  (+ (*field).u_comp[IND_2(n1i,k)]
					             - (*field).u_comp[IND_2(n1f-2,k)]
						     )*h_1;


	      }
	    }
	    

#ifdef EXCISION

	    /* interior and fases */
	  {register int n1;
	  register int n2;
	  for (n2 = exc_n2i; n2 < exc_n2f; n2++){
	      for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		(*dfield).u_comp[IND_2(n1,n2)] = 0.0;
	    }
	      (*dfield).u_comp[IND_2(exc_n1i-1,n2)] = dos*((*field).u_comp[IND_2(exc_n1i-1,n2)] 
		                                 - (*field).u_comp[IND_2(exc_n1i-2,n2)])*h_1;
	      (*dfield).u_comp[IND_2(exc_n1f,n2)] = dos*((*field).u_comp[IND_2(exc_n1f+1,n2)] 
		                                   - (*field).u_comp[IND_2(exc_n1f,n2)])*h_1;
	    }
	  }
	  

	    /* edges different from the one taking derivatives */

		(*dfield).u_comp[IND_2(exc_n1i-1,exc_n2i-1)] = - dos*(
                                 - (*field).u_comp[IND_2(exc_n1i+0,exc_n2i-1)] 
				 - (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-1)]        
			     + dos*(*field).u_comp[IND_2(exc_n1i-2,exc_n2i-1)]
				 )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1i-1,exc_n2f)] =  - dos*(
                                 - (*field).u_comp[IND_2(exc_n1i+0,exc_n2f)] 
				 - (*field).u_comp[IND_2(exc_n1i-1,exc_n2f)]        
			     + dos*(*field).u_comp[IND_2(exc_n1i-2,exc_n2f)]
				 )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2i-1)] =    dos*(
                               dos*(*field).u_comp[IND_2(exc_n1f+1,exc_n2i-1)] 
		                 - (*field).u_comp[IND_2(exc_n1f+0,exc_n2i-1)]        
				 - (*field).u_comp[IND_2(exc_n1f-1,exc_n2i-1)]
			       )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2f)] =    dos*(
                               dos*(*field).u_comp[IND_2(exc_n1f+1,exc_n2f)] 
			         - (*field).u_comp[IND_2(exc_n1f,exc_n2f)]        
			         - (*field).u_comp[IND_2(exc_n1f-1,exc_n2f)]
			       )*h_1/3.0;
	    


#endif	 



#ifdef DEBUG
	    printf("at case 1");
#endif
	    }

	break;


	/* Take derivative of second coordinate */

	case 2:{/* 	case 2: */


	    {register int k;
	    for (k = n1i; k < n1f; k++){
		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = dos*(- (*field).u_comp[IND_2(k,n2i)]
					           + (*field).u_comp[IND_2(k,n2i+1)]
						   )*h_2;



		/* Intermediate Values 1---(n-2) */

		{register int j;
		for (j = n2i+1; j < n2f-1; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (- (*field).u_comp[IND_2(k,j-1)]
					      	 + (*field).u_comp[IND_2(k,j+1)]
						 )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = dos*(+ (*field).u_comp[IND_2(k,n2f-1)]
					       	     - (*field).u_comp[IND_2(k,n2f-2)]
						     )*h_2;
	      }
	    }
	    

#ifdef EXCISION

	  {register int n1;
	  register int n2;
	    for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		for (n2 = exc_n2i; n2 < exc_n2f; n2++){
		  (*dfield).u_comp[IND_2(n1,n2)] = 0.0;
		}
		(*dfield).u_comp[IND_2(n1,exc_n2i-1)] = dos*((*field).u_comp[IND_2(n1,exc_n2i-1)] 
		                                   - (*field).u_comp[IND_2(n1,exc_n2i-2)])*h_2;
		(*dfield).u_comp[IND_2(n1,exc_n2f)] = dos*((*field).u_comp[IND_2(n1,exc_n2f+1)] 
						     - (*field).u_comp[IND_2(n1,exc_n2f)])*h_2;
	      }
	    }
	  

	    /* edges different from the one taking derivatives */



		(*dfield).u_comp[IND_2(exc_n1i-1,exc_n2i-1)] =  - dos*(
                                  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-0)] 
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-1)]        
				  + dos*(*field).u_comp[IND_2(exc_n1i-1,exc_n2i-2)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1i-1,exc_n2f)] =    dos*(
                                  + dos*(*field).u_comp[IND_2(exc_n1i-1,exc_n2f+1)] 
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2f+0)]        
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2f-1)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2i-1)] =  - dos*(
                                  -     (*field).u_comp[IND_2(exc_n1f,exc_n2i-0)] 
				  -     (*field).u_comp[IND_2(exc_n1f,exc_n2i-1)]        
				  + dos*(*field).u_comp[IND_2(exc_n1f,exc_n2i-2)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2f)] =    dos*(
                                    dos*(*field).u_comp[IND_2(exc_n1f,exc_n2f+1)] 
				    -  (*field).u_comp[IND_2(exc_n1f,exc_n2f)]        
				    -  (*field).u_comp[IND_2(exc_n1f,exc_n2f-1)]
				    )*h_2/3.0;
	    

  
#endif	 


#ifdef DEBUG
	    printf("at case 2");
#endif
	}
	
	break;


	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}

	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}

}


/* ------------------------------------------------------------------------------------------ */





/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


void derivS_Per_y_2d(struct grid_2d *grid,
	     struct field *field,
	     struct field *dfield,
	     int *coordinate)
{/* derivS_2d */

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;

    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;

#ifdef EXCISION
    int exc_n1 = (*grid).exc_grid_1;
    int exc_n2 = (*grid).exc_grid_2;
#endif

    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;


    FLOAT dos = 2.0;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi)*0.5;
    FLOAT h_2 = (FLOAT)(n2f-n2i)/(yf-yi)*0.5;     // Periodic

#ifdef EXCISION
    int exc_n1i = (n1f-n1i-1)/2-exc_n1+1; /* take out a -1 for visual compatibility with the rest of code */
    int exc_n2i = (n2f-n2i-1)/2-exc_n2+1;
    int exc_n1f = (n1f-n1i-1)/2+exc_n1; 
    int exc_n2f = (n2f-n2i-1)/2+exc_n2;
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{   /*  case 1: */

	    {register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = dos*(-(*field).u_comp[IND_2(n1i,k)]
					           +(*field).u_comp[IND_2(n1i+1,k)]
						   )*h_1;


		/* Intermediate Values 1---(n-2) */

		{register int j;
		for (j = n1i+1; j < n1f-1; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (- (*field).u_comp[IND_2(j-1,k)]
					         + (*field).u_comp[IND_2(j+1,k)]
						 )*h_1;
		}
		}

		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] = dos*(+ (*field).u_comp[IND_2(n1f-1,k)]
					             - (*field).u_comp[IND_2(n1f-2,k)]
						     )*h_1;


	      }
	    }
	    

#ifdef EXCISION

	    /* interior and fases */
	  {register int n1;
	  register int n2;
	  for (n2 = exc_n2i; n2 < exc_n2f; n2++){
	      for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		(*dfield).u_comp[IND_2(n1,n2)] = 0.0;
	    }
	      (*dfield).u_comp[IND_2(exc_n1i-1,n2)] = dos*((*field).u_comp[IND_2(exc_n1i-1,n2)] 
		                                 - (*field).u_comp[IND_2(exc_n1i-2,n2)])*h_1;
	      (*dfield).u_comp[IND_2(exc_n1f,n2)] = dos*((*field).u_comp[IND_2(exc_n1f+1,n2)] 
		                                   - (*field).u_comp[IND_2(exc_n1f,n2)])*h_1;
	    }
	  }
	  

	    /* edges different from the one taking derivatives */

		(*dfield).u_comp[IND_2(exc_n1i-1,exc_n2i-1)] = - dos*(
                                 - (*field).u_comp[IND_2(exc_n1i+0,exc_n2i-1)] 
				 - (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-1)]        
			     + dos*(*field).u_comp[IND_2(exc_n1i-2,exc_n2i-1)]
				 )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1i-1,exc_n2f)] =  - dos*(
                                 - (*field).u_comp[IND_2(exc_n1i+0,exc_n2f)] 
				 - (*field).u_comp[IND_2(exc_n1i-1,exc_n2f)]        
			     + dos*(*field).u_comp[IND_2(exc_n1i-2,exc_n2f)]
				 )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2i-1)] =    dos*(
                               dos*(*field).u_comp[IND_2(exc_n1f+1,exc_n2i-1)] 
		                 - (*field).u_comp[IND_2(exc_n1f+0,exc_n2i-1)]        
				 - (*field).u_comp[IND_2(exc_n1f-1,exc_n2i-1)]
			       )*h_1/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2f)] =    dos*(
                               dos*(*field).u_comp[IND_2(exc_n1f+1,exc_n2f)] 
			         - (*field).u_comp[IND_2(exc_n1f,exc_n2f)]        
			         - (*field).u_comp[IND_2(exc_n1f-1,exc_n2f)]
			       )*h_1/3.0;
	    


#endif	 



#ifdef DEBUG
	    printf("at case 1");
#endif
	    }

	break;


	/* Take derivative of second coordinate */

	case 2:{/* 	case 2: */


	    {register int k;
	    for (k = n1i; k < n1f; k++){
		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = (- (*field).u_comp[IND_2(k,n2f-1)]
					           + (*field).u_comp[IND_2(k,n2i+1)]
						   )*h_2;



		/* Intermediate Values 1---(n-2) */

		{register int j;
		for (j = n2i+1; j < n2f-1; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (- (*field).u_comp[IND_2(k,j-1)]
					      	 + (*field).u_comp[IND_2(k,j+1)]
						 )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = (+ (*field).u_comp[IND_2(k,n2i)]
					       	     - (*field).u_comp[IND_2(k,n2f-2)]
						     )*h_2;
	      }
	    }
	    

#ifdef EXCISION

	  {register int n1;
	  register int n2;
	    for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		for (n2 = exc_n2i; n2 < exc_n2f; n2++){
		  (*dfield).u_comp[IND_2(n1,n2)] = 0.0;
		}
		(*dfield).u_comp[IND_2(n1,exc_n2i-1)] = dos*((*field).u_comp[IND_2(n1,exc_n2i-1)] 
		                                   - (*field).u_comp[IND_2(n1,exc_n2i-2)])*h_2;
		(*dfield).u_comp[IND_2(n1,exc_n2f)] = dos*((*field).u_comp[IND_2(n1,exc_n2f+1)] 
						     - (*field).u_comp[IND_2(n1,exc_n2f)])*h_2;
	      }
	    }
	  

	    /* edges different from the one taking derivatives */



		(*dfield).u_comp[IND_2(exc_n1i-1,exc_n2i-1)] =  - dos*(
                                  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-0)] 
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2i-1)]        
				  + dos*(*field).u_comp[IND_2(exc_n1i-1,exc_n2i-2)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1i-1,exc_n2f)] =    dos*(
                                  + dos*(*field).u_comp[IND_2(exc_n1i-1,exc_n2f+1)] 
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2f+0)]        
				  -     (*field).u_comp[IND_2(exc_n1i-1,exc_n2f-1)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2i-1)] =  - dos*(
                                  -     (*field).u_comp[IND_2(exc_n1f,exc_n2i-0)] 
				  -     (*field).u_comp[IND_2(exc_n1f,exc_n2i-1)]        
				  + dos*(*field).u_comp[IND_2(exc_n1f,exc_n2i-2)]
				  )*h_2/3.0;

                (*dfield).u_comp[IND_2(exc_n1f,exc_n2f)] =    dos*(
                                    dos*(*field).u_comp[IND_2(exc_n1f,exc_n2f+1)] 
				    -  (*field).u_comp[IND_2(exc_n1f,exc_n2f)]        
				    -  (*field).u_comp[IND_2(exc_n1f,exc_n2f-1)]
				    )*h_2/3.0;
	    

  
#endif	 


#ifdef DEBUG
	    printf("at case 2");
#endif
	}
	
	break;


	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}

	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}

}





/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */







void derivQ_2d(struct grid_2d *grid,
	     struct field *field,
	     struct field *dfield,
	     int *coordinate)
{ /* derivQ_2d */

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;
    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;

#ifdef EXCISION

    int exc_n1 = (*grid).exc_grid_1;
    int exc_n2 = (*grid).exc_grid_2;
    int exc_n1i = (n1f-n1i-1)/2-exc_n1+1; /* a +1 for visual compatibility with the rest of code */
    int exc_n2i = (n2f-n2i-1)/2-exc_n2+1;
    int exc_n1f = (n1f-n1i-1)/2+exc_n1; 
    int exc_n2f = (n2f-n2i-1)/2+exc_n2;
#endif

    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;

    FLOAT d00 = -1.4117647059;
    FLOAT d01 =  1.7352941176;
    FLOAT d02 = -0.23529411765;
    FLOAT d03 = -0.088235294118;

    FLOAT d10 = -0.5;
    FLOAT d12 = -d10;

    FLOAT d20 = 0.093023255814;
    FLOAT d21 = -0.68604651163;
    FLOAT d23 = -d21;
    FLOAT d24 = -d20;

    FLOAT d30 = 0.030612244898;
    FLOAT d32 = -0.60204081633;
    FLOAT d34 = 0.65306122449;
    FLOAT d35 = -0.081632653061;

    FLOAT d1 = 0.6666666666666;
    FLOAT d2 = -0.0833333333333;


    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);

#ifdef EXCISION

    printf("Warning with excision this is not correct at corners and edges");
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{    /* case 1: */
	    {register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = (d00*(*field).u_comp[IND_2(n1i,k)] 
					    + d01*(*field).u_comp[IND_2(n1i+1,k)]
					    + d02*(*field).u_comp[IND_2(n1i+2,k)] 
					    + d03*(*field).u_comp[IND_2(n1i+3,k)]
					    )*h_1;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(n1i+1,k)] = (d10*(*field).u_comp[IND_2(n1i+0,k)]
					      + d12*(*field).u_comp[IND_2(n1i+2,k)]
					      )*h_1;

		/* Value at 2 */

		(*dfield).u_comp[IND_2(n1i+2,k)] = (d20*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d21*(*field).u_comp[IND_2(n1i+1,k)]
					      + d23*(*field).u_comp[IND_2(n1i+3,k)] 
					      + d24*(*field).u_comp[IND_2(n1i+4,k)]
					      )*h_1;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(n1i+3,k)] = (d30*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d32*(*field).u_comp[IND_2(n1i+2,k)]
					      + d34*(*field).u_comp[IND_2(n1i+4,k)] 
					      + d35*(*field).u_comp[IND_2(n1i+5,k)]
					      )*h_1;

		/* Intermediate Values 4---(n-5) */

		{register int j; 
		for (j = n1i+4; j < n1f-4; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (-d2*(*field).u_comp[IND_2(j-2,k)] 
					      - d1*(*field).u_comp[IND_2(j-1,k)]
					      + d1*(*field).u_comp[IND_2(j+1,k)] 
					      + d2*(*field).u_comp[IND_2(j+2,k)]
					      )*h_1;
		}
		}

		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] = -(d00*(*field).u_comp[IND_2(n1f-1,k)] 
					       + d01*(*field).u_comp[IND_2(n1f-2,k)]
					       + d02*(*field).u_comp[IND_2(n1f-3,k)]
					       + d03*(*field).u_comp[IND_2(n1f-4,k)]
					       )*h_1;

		/* Value at n-2 */

		(*dfield).u_comp[IND_2(n1f-2,k)] = -(d10*(*field).u_comp[IND_2(n1f-1,k)]
					       + d12*(*field).u_comp[IND_2(n1f-3,k)]
					       )*h_1;

		/* Value at n-3 */

		(*dfield).u_comp[IND_2(n1f-3,k)] = -( d20*(*field).u_comp[IND_2(n1f-1,k)]
						+ d21*(*field).u_comp[IND_2(n1f-2,k)]
						+ d23*(*field).u_comp[IND_2(n1f-4,k)]
						+ d24*(*field).u_comp[IND_2(n1f-5,k)]
						)*h_1;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(n1f-4,k)] = -(d30*(*field).u_comp[IND_2(n1f-1,k)]
					       + d32*(*field).u_comp[IND_2(n1f-3,k)]
					       + d34*(*field).u_comp[IND_2(n1f-5,k)]
					       + d35*(*field).u_comp[IND_2(n1f-6,k)]
					       )*h_1;
	      }
	    }
	    


#ifdef EXCISION
#ifdef DEBUG
    printf("exc_n1 = %d, exc_n1i = %d, exc_n1f = %d \n",   exc_n1, exc_n1i, exc_n1f);
#endif
	  {register int n1;
	  register int n2;
	  for (n2 = exc_n2i; n2 < exc_n2f; n2++){
	      for (n1 = exc_n1i; n1 < exc_n1f; n1++){
		(*dfield).u_comp[IND_2(n1,n2)] = 0.0;
#ifdef DEBUG
    printf("n1 = %d, n2 =  %d \n",   n1, n2);
#endif
	    }
	  
	  



	        (*dfield).u_comp[IND_2(exc_n1f,n2)] = (  d00*(*field).u_comp[IND_2(exc_n1f,n2)] 
						     + d01*(*field).u_comp[IND_2(exc_n1f+1,n2)]
						     + d02*(*field).u_comp[IND_2(exc_n1f+2,n2)] 
						     + d03*(*field).u_comp[IND_2(exc_n1f+3,n2)]
						     )*h_1;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(exc_n1f+1,n2)] = (  d10*(*field).u_comp[IND_2(exc_n1f+0,n2)]
						       + d12*(*field).u_comp[IND_2(exc_n1f+2,n2)]
						       )*h_1;

		/* Value at 2 */

		(*dfield).u_comp[IND_2(exc_n1f+2,n2)] = (  d20*(*field).u_comp[IND_2(exc_n1f+0,n2)] 
						       + d21*(*field).u_comp[IND_2(exc_n1f+1,n2)]
						       + d23*(*field).u_comp[IND_2(exc_n1f+3,n2)] 
						       + d24*(*field).u_comp[IND_2(exc_n1f+4,n2)]
						       )*h_1;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(exc_n1f+3,n2)] = (  d30*(*field).u_comp[IND_2(exc_n1f+0,n2)] 
						       + d32*(*field).u_comp[IND_2(exc_n1f+2,n2)]
						       + d34*(*field).u_comp[IND_2(exc_n1f+4,n2)] 
						       + d35*(*field).u_comp[IND_2(exc_n1f+5,n2)]
						       )*h_1;
		/* Value at n-1 */


		(*dfield).u_comp[IND_2(exc_n1i-1,n2)] = -(  d00*(*field).u_comp[IND_2(exc_n1i-1,n2)] 
							+ d01*(*field).u_comp[IND_2(exc_n1i-2,n2)]
							+ d02*(*field).u_comp[IND_2(exc_n1i-3,n2)]
							+ d03*(*field).u_comp[IND_2(exc_n1i-4,n2)]
							)*h_1;

		/* Value at n-2 */

		(*dfield).u_comp[IND_2(exc_n1i-2,n2)] = -(  d10*(*field).u_comp[IND_2(exc_n1i-1,n2)]
							+ d12*(*field).u_comp[IND_2(exc_n1i-3,n2)]
							)*h_1;

		/* Value at n-3 */

		(*dfield).u_comp[IND_2(exc_n1i-3,n2)] = -(   d20*(*field).u_comp[IND_2(exc_n1i-1,n2)]
							 + d21*(*field).u_comp[IND_2(exc_n1i-2,n2)]
							 + d23*(*field).u_comp[IND_2(exc_n1i-4,n2)]
							 + d24*(*field).u_comp[IND_2(exc_n1i-5,n2)]
							 )*h_1;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(exc_n1i-4,n2)] = -(  d30*(*field).u_comp[IND_2(exc_n1i-1,n2)]
							+ d32*(*field).u_comp[IND_2(exc_n1i-3,n2)]
							+ d34*(*field).u_comp[IND_2(exc_n1i-5,n2)]
							+ d35*(*field).u_comp[IND_2(exc_n1i-6,n2)]
							)*h_1;
	    }
	  }
	  
#endif	 



#ifdef DEBUG
	    printf("at case 1");
#endif
	}
	    
	
	break;
	

	/* Take derivative of second coordinate */

	case 2:{   /* case 2: */
	    
	    {register int k;
	    register int l;
	    for (k = n1i; k < n1f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = (d00*(*field).u_comp[IND_2(k,n2i)]
					    + d01*(*field).u_comp[IND_2(k,n2i+1)]
					    + d02*(*field).u_comp[IND_2(k,n2i+2)]
					    + d03*(*field).u_comp[IND_2(k,n2i+3)]
					    )*h_2;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(k,n2i+1)] = (d10*(*field).u_comp[IND_2(k,n2i+0)]
					      + d12*(*field).u_comp[IND_2(k,n2i+2)]
					      )*h_2;
		
		/* Value at 2 */
    
		(*dfield).u_comp[IND_2(k,n2i+2)] = (d20*(*field).u_comp[IND_2(k,n2i+0)]
					      + d21*(*field).u_comp[IND_2(k,n2i+1)]
					      + d23*(*field).u_comp[IND_2(k,n2i+3)]
					      + d24*(*field).u_comp[IND_2(k,n2i+4)]
					      )*h_2;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(k,n2i+3)] = (d30*(*field).u_comp[IND_2(k,n2i+0)]
					      + d32*(*field).u_comp[IND_2(k,n2i+2)]
					      + d34*(*field).u_comp[IND_2(k,n2i+4)]
					      + d35*(*field).u_comp[IND_2(k,n2i+5)]
					      )*h_2;

		/* Intermediate Values 4---(n-5) */

		{register int j;
		for (j = n2i+4; j < n2f-4; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (-d2*(*field).u_comp[IND_2(k,j-2)]
					      - d1*(*field).u_comp[IND_2(k,j-1)]
					      + d1*(*field).u_comp[IND_2(k,j+1)]
					      + d2*(*field).u_comp[IND_2(k,j+2)]
					      )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = -(d00*(*field).u_comp[IND_2(k,n2f-1)]
					       + d01*(*field).u_comp[IND_2(k,n2f-2)]
					       + d02*(*field).u_comp[IND_2(k,n2f-3)]
					       + d03*(*field).u_comp[IND_2(k,n2f-4)]
					       )*h_2;
		
		/* Value at n-2 */
		
		(*dfield).u_comp[IND_2(k,n2f-2)] = -(d10*(*field).u_comp[IND_2(k,n2f-1)]
					       + d12*(*field).u_comp[IND_2(k,n2f-3)]
					       )*h_2;
		
		/* Value at n-3 */
		
		(*dfield).u_comp[IND_2(k,n2f-3)] = -( d20*(*field).u_comp[IND_2(k,n2f-1)]
						+ d21*(*field).u_comp[IND_2(k,n2f-2)]
						+ d23*(*field).u_comp[IND_2(k,n2f-4)]
						+ d24*(*field).u_comp[IND_2(k,n2f-5)]
						)*h_2;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(k,n2f-4)] = -(d30*(*field).u_comp[IND_2(k,n2f-1)]
					       + d32*(*field).u_comp[IND_2(k,n2f-3)]
					       + d34*(*field).u_comp[IND_2(k,n2f-5)]
					       + d35*(*field).u_comp[IND_2(k,n2f-6)]
					       )*h_2;
	      }
	    }
	    


#ifdef EXCISION       

	  {register int n1;
	  register int n2;
	  for (n1 = exc_n1i; n1 < exc_n1f; n1++){
	      for (n2 = exc_n2i; n2 < exc_n2f; n2++){
		(*dfield).u_comp[IND_2(n1,n2)] = 0.0;
	    }



	        (*dfield).u_comp[IND_2(n1,exc_n2f)] = (  d00*(*field).u_comp[IND_2(n1,exc_n2f)] 
						     + d01*(*field).u_comp[IND_2(n1,exc_n2f+1)]
						     + d02*(*field).u_comp[IND_2(n1,exc_n2f+2)] 
						     + d03*(*field).u_comp[IND_2(n1,exc_n2f+3)]
						     )*h_2;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(n1,exc_n2f+1)] = (  d10*(*field).u_comp[IND_2(n1,exc_n2f+0)]
						       + d12*(*field).u_comp[IND_2(n1,exc_n2f+2)]
						       )*h_2;

		/* Value at 2 */

		(*dfield).u_comp[IND_2(n1,exc_n2f+2)] = (  d20*(*field).u_comp[IND_2(n1,exc_n2f+0)] 
						       + d21*(*field).u_comp[IND_2(n1,exc_n2f+1)]
						       + d23*(*field).u_comp[IND_2(n1,exc_n2f+3)] 
						       + d24*(*field).u_comp[IND_2(n1,exc_n2f+4)]
						       )*h_2;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(n1,exc_n2f+3)] = (  d30*(*field).u_comp[IND_2(n1,exc_n2f+0)] 
						       + d32*(*field).u_comp[IND_2(n1,exc_n2f+2)]
						       + d34*(*field).u_comp[IND_2(n1,exc_n2f+4)] 
						       + d35*(*field).u_comp[IND_2(n1,exc_n2f+5)]
						       )*h_2;
		/* Value at n-1 */


		(*dfield).u_comp[IND_2(n1,exc_n2i-1)] = -(  d00*(*field).u_comp[IND_2(n1,exc_n2i-1)] 
							+ d01*(*field).u_comp[IND_2(n1,exc_n2i-2)]
							+ d02*(*field).u_comp[IND_2(n1,exc_n2i-3)]
							+ d03*(*field).u_comp[IND_2(n1,exc_n2i-4)]
							)*h_2;

		/* Value at n-2 */

		(*dfield).u_comp[IND_2(n1,exc_n2i-2)] = -(  d10*(*field).u_comp[IND_2(n1,exc_n2i-1)]
							+ d12*(*field).u_comp[IND_2(n1,exc_n2i-3)]
							)*h_2;

		/* Value at n-3 */

		(*dfield).u_comp[IND_2(n1,exc_n2i-3)] = -(   d20*(*field).u_comp[IND_2(n1,exc_n2i-1)]
							 + d21*(*field).u_comp[IND_2(n1,exc_n2i-2)]
							 + d23*(*field).u_comp[IND_2(n1,exc_n2i-4)]
							 + d24*(*field).u_comp[IND_2(n1,exc_n2i-5)]
							 )*h_2;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(n1,exc_n2i-4)] = -(  d30*(*field).u_comp[IND_2(n1,exc_n2i-1)]
							+ d32*(*field).u_comp[IND_2(n1,exc_n2i-3)]
							+ d34*(*field).u_comp[IND_2(n1,exc_n2i-5)]
							+ d35*(*field).u_comp[IND_2(n1,exc_n2i-6)]
							)*h_2;
	    }
	  }
	  
#endif	 



#ifdef DEBUG
	    printf("at case 2");
#endif
	}
	
	break;

	
	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}

	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}
}





/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */



/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */



void derivQ_Per_y_2d(struct grid_2d *grid,
	     struct field *field,
	     struct field *dfield,
	     int *coordinate)
{ /* derivQ_Per_y_2d */

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;
    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;
    
    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;

    FLOAT d00 = -1.4117647058823529411764705882353;
    FLOAT d01 = 1.7352941176470588235294117647059;
    FLOAT d02 = -0.23529411764705882352941176470588;
    FLOAT d03 = -0.088235294117647058823529411764706;

    FLOAT d10 = -0.50000000000000000000000000000000;
    FLOAT d12 = 0.50000000000000000000000000000000;

    FLOAT d20 = 0.093023255813953488372093023255814;
    FLOAT d21 = -0.68604651162790697674418604651163;
    FLOAT d23 = 0.68604651162790697674418604651163;
    FLOAT d24 = -0.093023255813953488372093023255814;

    FLOAT d30 = 0.030612244897959183673469387755102;
    FLOAT d32 = -0.60204081632653061224489795918367;
    FLOAT d34 = 0.65306122448979591836734693877551;
    FLOAT d35 = -0.081632653061224489795918367346939;

    FLOAT d1 = 0.6666666666666;
    FLOAT d2 = -0.0833333333333;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
    FLOAT h_2 = (FLOAT)(n2f-n2i)/(yf-yi); // for periodic removed the 1

#ifdef EXCISION
		        printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION ! **************");
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{    /* case 1: */
	    {register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = (d00*(*field).u_comp[IND_2(n1i,k)] 
					    + d01*(*field).u_comp[IND_2(n1i+1,k)]
					    + d02*(*field).u_comp[IND_2(n1i+2,k)] 
					    + d03*(*field).u_comp[IND_2(n1i+3,k)]
					    )*h_1;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(n1i+1,k)] = (d10*(*field).u_comp[IND_2(n1i+0,k)]
					      + d12*(*field).u_comp[IND_2(n1i+2,k)]
					      )*h_1;

		/* Value at 2 */

		(*dfield).u_comp[IND_2(n1i+2,k)] = (d20*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d21*(*field).u_comp[IND_2(n1i+1,k)]
					      + d23*(*field).u_comp[IND_2(n1i+3,k)] 
					      + d24*(*field).u_comp[IND_2(n1i+4,k)]
					      )*h_1;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(n1i+3,k)] = (d30*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d32*(*field).u_comp[IND_2(n1i+2,k)]
					      + d34*(*field).u_comp[IND_2(n1i+4,k)] 
					      + d35*(*field).u_comp[IND_2(n1i+5,k)]
					      )*h_1;

		/* Intermediate Values 4---(n-5) */

		{register int j; 
		for (j = n1i+4; j < n1f-4; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (-d2*(*field).u_comp[IND_2(j-2,k)] 
					      - d1*(*field).u_comp[IND_2(j-1,k)]
					      + d1*(*field).u_comp[IND_2(j+1,k)] 
					      + d2*(*field).u_comp[IND_2(j+2,k)]
					      )*h_1;
		}
		}

		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] = -(d00*(*field).u_comp[IND_2(n1f-1,k)] 
					       + d01*(*field).u_comp[IND_2(n1f-2,k)]
					       + d02*(*field).u_comp[IND_2(n1f-3,k)]
					       + d03*(*field).u_comp[IND_2(n1f-4,k)]
					       )*h_1;

		/* Value at n-2 */

		(*dfield).u_comp[IND_2(n1f-2,k)] = -(d10*(*field).u_comp[IND_2(n1f-1,k)]
					       + d12*(*field).u_comp[IND_2(n1f-3,k)]
					       )*h_1;

		/* Value at n-3 */

		(*dfield).u_comp[IND_2(n1f-3,k)] = -( d20*(*field).u_comp[IND_2(n1f-1,k)]
						+ d21*(*field).u_comp[IND_2(n1f-2,k)]
						+ d23*(*field).u_comp[IND_2(n1f-4,k)]
						+ d24*(*field).u_comp[IND_2(n1f-5,k)]
						)*h_1;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(n1f-4,k)] = -(d30*(*field).u_comp[IND_2(n1f-1,k)]
					       + d32*(*field).u_comp[IND_2(n1f-3,k)]
					       + d34*(*field).u_comp[IND_2(n1f-5,k)]
					       + d35*(*field).u_comp[IND_2(n1f-6,k)]
					       )*h_1;
	      }
	    }
	    



#ifdef DEBUG
	    printf("at case 1");
#endif
	}
	    
	
	break;
	

	/* Take derivative of second coordinate THIS IS PERIODIC */

	case 2:{   /* case 2: */
	    
	    {register int k;
	    register int l;
	    for (k = n1i; k < n1f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = (-d2*(*field).u_comp[IND_2(k,n2f-2)]
					          -d1*(*field).u_comp[IND_2(k,n2f-1)]
					          +d1*(*field).u_comp[IND_2(k,n2i+1)]
					          +d2*(*field).u_comp[IND_2(k,n2i+2)]
					         )*h_2;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(k,n2i+1)] = (-d2*(*field).u_comp[IND_2(k,n2f-1)]
					            -d1*(*field).u_comp[IND_2(k,n2i-0)]
					            +d1*(*field).u_comp[IND_2(k,n2i+2)]
					            +d2*(*field).u_comp[IND_2(k,n2i+3)]
					           )*h_2;
		

		/* Intermediate Values 2---(n-3) */

		{register int j;
		for (j = n2i+2; j < n2f-2; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (-d2*(*field).u_comp[IND_2(k,j-2)]
					      - d1*(*field).u_comp[IND_2(k,j-1)]
					      + d1*(*field).u_comp[IND_2(k,j+1)]
					      + d2*(*field).u_comp[IND_2(k,j+2)]
					      )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = (-d2*(*field).u_comp[IND_2(k,n2f-3)]
					            -d1*(*field).u_comp[IND_2(k,n2f-2)]
					            +d1*(*field).u_comp[IND_2(k,n2i+0)]
					            +d2*(*field).u_comp[IND_2(k,n2i+1)]
					           )*h_2;
		
		
		/* Value at n-2 */
		
		(*dfield).u_comp[IND_2(k,n2f-2)] = (-d2*(*field).u_comp[IND_2(k,n2f-4)]
					            -d1*(*field).u_comp[IND_2(k,n2f-3)]
					            +d1*(*field).u_comp[IND_2(k,n2f-1)]
					            +d2*(*field).u_comp[IND_2(k,n2i+0)]
					           )*h_2;
	      }
	    }


#ifdef DEBUG
	    printf("at case 2");
#endif
	}
	
	break;

	
	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}

	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}
}

/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */



void derivQ_Per_x_2d(struct grid_2d *grid,
	     struct field *field,
	     struct field *dfield,
	     int *coordinate)
{ // derivQ_Per_x_2d 

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;
    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;
    
    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;

    FLOAT d00 = -1.4117647058823529411764705882353;
    FLOAT d01 = 1.7352941176470588235294117647059;
    FLOAT d02 = -0.23529411764705882352941176470588;
    FLOAT d03 = -0.088235294117647058823529411764706;

    FLOAT d10 = -0.50000000000000000000000000000000;
    FLOAT d12 = 0.50000000000000000000000000000000;

    FLOAT d20 = 0.093023255813953488372093023255814;
    FLOAT d21 = -0.68604651162790697674418604651163;
    FLOAT d23 = 0.68604651162790697674418604651163;
    FLOAT d24 = -0.093023255813953488372093023255814;

    FLOAT d30 = 0.030612244897959183673469387755102;
    FLOAT d32 = -0.60204081632653061224489795918367;
    FLOAT d34 = 0.65306122448979591836734693877551;
    FLOAT d35 = -0.081632653061224489795918367346939;
    
    FLOAT d1 = 0.6666666666666;
    FLOAT d2 = -0.0833333333333;

    FLOAT h_1 = (FLOAT)(n1f-n1i)/(xf-xi);         // for periodic removed the 1
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi); 

#ifdef EXCISION
		        printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION ! **************");
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{    /* case 1: */
	    {register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = (-d2*(*field).u_comp[IND_2(n1f-2,k)] 
					          -d1*(*field).u_comp[IND_2(n1f-1,k)]
					          +d1*(*field).u_comp[IND_2(n1i+1,k)] 
					          +d2*(*field).u_comp[IND_2(n1i+2,k)]
					         )*h_1;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(n1i+1,k)] = (-d2*(*field).u_comp[IND_2(n1f-1,k)] 
					            -d1*(*field).u_comp[IND_2(n1i+0,k)]
					            +d1*(*field).u_comp[IND_2(n1i+2,k)] 
					            +d2*(*field).u_comp[IND_2(n1i+3,k)]
					            )*h_1;


		/* Intermediate Values 2---(n-3) */

		{register int j; 
		for (j = n1i+2; j < n1f-2; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (-d2*(*field).u_comp[IND_2(j-2,k)] 
					            -d1*(*field).u_comp[IND_2(j-1,k)]
					            +d1*(*field).u_comp[IND_2(j+1,k)] 
					            +d2*(*field).u_comp[IND_2(j+2,k)]
					            )*h_1;
		}
		}

		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] = (-d2*(*field).u_comp[IND_2(n1f-3,k)] 
					            -d1*(*field).u_comp[IND_2(n1f-2,k)]
					            +d1*(*field).u_comp[IND_2(n1i+0,k)] 
					            +d2*(*field).u_comp[IND_2(n1i+1,k)]
					            )*h_1;


		/* Value at n-2 */

		(*dfield).u_comp[IND_2(n1f-2,k)] = (-d2*(*field).u_comp[IND_2(n1f-4,k)] 
					            -d1*(*field).u_comp[IND_2(n1f-3,k)]
					            +d1*(*field).u_comp[IND_2(n1f-1,k)] 
					            +d2*(*field).u_comp[IND_2(n1i+0,k)]
					            )*h_1;

	
	      }
	    }
	    



#ifdef DEBUG
	    printf("at case 1");
#endif
	}
	    
	
	break;
	
	/* Take derivative of second coordinate */
	


	case 2:{   /* case 2: */
	    
	    {register int k;
	    register int l;
	    for (k = n1i; k < n1f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = (d00*(*field).u_comp[IND_2(k,n2i)]
					    + d01*(*field).u_comp[IND_2(k,n2i+1)]
					    + d02*(*field).u_comp[IND_2(k,n2i+2)]
					    + d03*(*field).u_comp[IND_2(k,n2i+3)]
					    )*h_2;

		/* Value at 1 */

		(*dfield).u_comp[IND_2(k,n2i+1)] = (d10*(*field).u_comp[IND_2(k,n2i+0)]
					      + d12*(*field).u_comp[IND_2(k,n2i+2)]
					      )*h_2;
		
		/* Value at 2 */
    
		(*dfield).u_comp[IND_2(k,n2i+2)] = (d20*(*field).u_comp[IND_2(k,n2i+0)]
					      + d21*(*field).u_comp[IND_2(k,n2i+1)]
					      + d23*(*field).u_comp[IND_2(k,n2i+3)]
					      + d24*(*field).u_comp[IND_2(k,n2i+4)]
					      )*h_2;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(k,n2i+3)] = (d30*(*field).u_comp[IND_2(k,n2i+0)]
					      + d32*(*field).u_comp[IND_2(k,n2i+2)]
					      + d34*(*field).u_comp[IND_2(k,n2i+4)]
					      + d35*(*field).u_comp[IND_2(k,n2i+5)]
					      )*h_2;

		/* Intermediate Values 4---(n-5) */

		{register int j;
		for (j = n2i+4; j < n2f-4; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (-d2*(*field).u_comp[IND_2(k,j-2)]
					      - d1*(*field).u_comp[IND_2(k,j-1)]
					      + d1*(*field).u_comp[IND_2(k,j+1)]
					      + d2*(*field).u_comp[IND_2(k,j+2)]
					      )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = -(d00*(*field).u_comp[IND_2(k,n2f-1)]
					       + d01*(*field).u_comp[IND_2(k,n2f-2)]
					       + d02*(*field).u_comp[IND_2(k,n2f-3)]
					       + d03*(*field).u_comp[IND_2(k,n2f-4)]
					       )*h_2;
		
		/* Value at n-2 */
		
		(*dfield).u_comp[IND_2(k,n2f-2)] = -(d10*(*field).u_comp[IND_2(k,n2f-1)]
					       + d12*(*field).u_comp[IND_2(k,n2f-3)]
					       )*h_2;
		
		/* Value at n-3 */
		
		(*dfield).u_comp[IND_2(k,n2f-3)] = -( d20*(*field).u_comp[IND_2(k,n2f-1)]
						+ d21*(*field).u_comp[IND_2(k,n2f-2)]
						+ d23*(*field).u_comp[IND_2(k,n2f-4)]
						+ d24*(*field).u_comp[IND_2(k,n2f-5)]
						)*h_2;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(k,n2f-4)] = -(d30*(*field).u_comp[IND_2(k,n2f-1)]
					       + d32*(*field).u_comp[IND_2(k,n2f-3)]
					       + d34*(*field).u_comp[IND_2(k,n2f-5)]
					       + d35*(*field).u_comp[IND_2(k,n2f-6)]
					       )*h_2;
	      }
	    }
	    



#ifdef DEBUG
	    printf("at case 2");
#endif
	}
	
	break;

	
	
	
	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}

	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}
}

/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */







void derivQ_3_2d(struct grid_2d *grid,
	     struct field *field, 
	     struct field *dfield, 
	     int *coordinate)
{   /* derivQ_3_2d */

    int n1 = (*grid).n_gridpts_1;
    int n2 = (*grid).n_gridpts_2;
    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1;
    int n_gridpts_2 = n2;

    int coord = (*coordinate);

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;

    FLOAT d00 = -1.833333333333333333333333;
    FLOAT d01 =  3.0;
    FLOAT d02 =  -1.5;
    FLOAT d03 =  0.333333333333333333333333;
    
    FLOAT d10 = -0.38942207148531184298;
    FLOAT d11 = -0.26953763903486946056;
    FLOAT d12 =  0.63903793765926293838;
    FLOAT d13 =  0.094332736084546377480;
    FLOAT d14 = -0.080518371580844513359;
    FLOAT d15 =  0.0061074083572165009295;

    FLOAT d20 =  0.11124996667625322721;
    FLOAT d21 = -0.78615310943278550936;
    FLOAT d22 =  0.19877943763527643222;
    FLOAT d23 =  0.50808067692835148792;
    FLOAT d24 = -0.024137062412656370601;
    FLOAT d25 = -0.0078199093944392672116;

    FLOAT d30 =  0.019051206094885019047822;
    FLOAT d31 =  0.026931104200732614181666;
    FLOAT d32 = -0.633860292039252305642283;
    FLOAT d33 =  0.051772670918649366462688;
    FLOAT d34 =  0.592764606048964306931634;
    FLOAT d35 = -0.054368814269840675877468;
    FLOAT d36 = -0.002290480954138351040607;

    FLOAT d40 = -0.002498706495423627386248;
    FLOAT d41 =  0.005463924453044550084942;
    FLOAT d42 =  0.087024805619019315445041;
    FLOAT d43 = -0.686097670431383548237962; 
    FLOAT d44 =  0.018985530480943661987934;
    FLOAT d45 =  0.659895344563505072850627;
    FLOAT d46 = -0.082773228189705424744336;


    FLOAT d1 = 0.6666666666666666666666666;
    FLOAT d2 = -0.0833333333333333333333333;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);

#ifdef EXCISION
    printf("Waning excision not implemented in derivQ_3_2d");
#endif

    switch(coord)
	{

   /* Take derivative of first coordinate */

	case 1:{   /* case 1: */
	    register int k;
	    for (k = n2i; k < n2f; k++){

		/* Value at 0 */

		(*dfield).u_comp[IND_2(n1i,k)] = (d00*(*field).u_comp[IND_2(n1i,k)] 
					    + d01*(*field).u_comp[IND_2(n1i+1,k)]
					    + d02*(*field).u_comp[IND_2(n1i+2,k)] 
					    + d03*(*field).u_comp[IND_2(n1i+3,k)]
					    )*h_1;
		
		/* Value at 1 */
		
		(*dfield).u_comp[IND_2(n1i+1,k)] = (d10*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d11*(*field).u_comp[IND_2(n1i+1,k)]
					      + d12*(*field).u_comp[IND_2(n1i+2,k)]
					      + d13*(*field).u_comp[IND_2(n1i+3,k)]
					      + d14*(*field).u_comp[IND_2(n1i+4,k)]
					      + d15*(*field).u_comp[IND_2(n1i+5,k)]
					      )*h_1;
  	
		/* Value at 2 */
    
		(*dfield).u_comp[IND_2(n1i+2,k)] = (d20*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d21*(*field).u_comp[IND_2(n1i+1,k)]
					      + d22*(*field).u_comp[IND_2(n1i+2,k)]
					      + d23*(*field).u_comp[IND_2(n1i+3,k)] 
					      + d24*(*field).u_comp[IND_2(n1i+4,k)]
					      + d25*(*field).u_comp[IND_2(n1i+5,k)]
					      )*h_1;

		/* Value at 3 */
		
		(*dfield).u_comp[IND_2(n1i+3,k)] = (d30*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d31*(*field).u_comp[IND_2(n1i+1,k)]
					      + d32*(*field).u_comp[IND_2(n1i+2,k)]
					      + d33*(*field).u_comp[IND_2(n1i+3,k)]
					      + d34*(*field).u_comp[IND_2(n1i+4,k)] 
					      + d35*(*field).u_comp[IND_2(n1i+5,k)]
					      + d36*(*field).u_comp[IND_2(n1i+6,k)]
					      )*h_1;

		/* Value at 4 */
		
		(*dfield).u_comp[IND_2(n1i+4,k)] = (d40*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d41*(*field).u_comp[IND_2(n1i+1,k)]
					      + d42*(*field).u_comp[IND_2(n1i+2,k)]
					      + d43*(*field).u_comp[IND_2(n1i+3,k)]
					      + d44*(*field).u_comp[IND_2(n1i+4,k)] 
					      + d45*(*field).u_comp[IND_2(n1i+5,k)]
					      + d46*(*field).u_comp[IND_2(n1i+6,k)]
					      )*h_1;

		/* Intermediate Values 5---(n-6) */

		{register int j; 
		for (j = n1i+5; j < n1f-5; ++j){
		    (*dfield).u_comp[IND_2(j,k)] = (-d2*(*field).u_comp[IND_2(j-2,k)] 
					      - d1*(*field).u_comp[IND_2(j-1,k)]
					      + d1*(*field).u_comp[IND_2(j+1,k)] 
					      + d2*(*field).u_comp[IND_2(j+2,k)]
					      )*h_1;
		}
		}
 
		/* Value at n-1 */

		(*dfield).u_comp[IND_2(n1f-1,k)] = -(d00*(*field).u_comp[IND_2(n1f-1,k)] 
					       + d01*(*field).u_comp[IND_2(n1f-2,k)]
					       + d02*(*field).u_comp[IND_2(n1f-3,k)] 
					       + d03*(*field).u_comp[IND_2(n1f-4,k)]
					       )*h_1;
    
		/* Value at n-2 */
		
		(*dfield).u_comp[IND_2(n1f-2,k)] = -(d10*(*field).u_comp[IND_2(n1f-1,k)] 
					       + d11*(*field).u_comp[IND_2(n1f-2,k)]
					       + d12*(*field).u_comp[IND_2(n1f-3,k)]
					       + d13*(*field).u_comp[IND_2(n1f-4,k)]
					       + d14*(*field).u_comp[IND_2(n1f-5,k)]
					       + d15*(*field).u_comp[IND_2(n1f-6,k)]
					       )*h_1;
		
		/* Value at n-3 */
    
		(*dfield).u_comp[IND_2(n1f-3,k)] = -( d20*(*field).u_comp[IND_2(n1f-1,k)] 
						+ d21*(*field).u_comp[IND_2(n1f-2,k)]
						+ d22*(*field).u_comp[IND_2(n1f-3,k)]
						+ d23*(*field).u_comp[IND_2(n1f-4,k)]
						+ d24*(*field).u_comp[IND_2(n1f-5,k)]
						+ d25*(*field).u_comp[IND_2(n1f-6,k)]
						)*h_1;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(n1f-4,k)] = -(d30*(*field).u_comp[IND_2(n1f-1,k)]
					       + d31*(*field).u_comp[IND_2(n1f-2,k)]
					       + d32*(*field).u_comp[IND_2(n1f-3,k)]
					       + d33*(*field).u_comp[IND_2(n1f-4,k)]
					       + d34*(*field).u_comp[IND_2(n1f-5,k)] 
					       + d35*(*field).u_comp[IND_2(n1f-6,k)]
					       + d36*(*field).u_comp[IND_2(n1f-7,k)]
					       )*h_1;
	    	/* Value at n-5 */

		(*dfield).u_comp[IND_2(n1f-5,k)] = -(d40*(*field).u_comp[IND_2(n1f-1,k)]
					       + d41*(*field).u_comp[IND_2(n1f-2,k)]
					       + d42*(*field).u_comp[IND_2(n1f-3,k)]
					       + d43*(*field).u_comp[IND_2(n1f-4,k)]
					       + d44*(*field).u_comp[IND_2(n1f-5,k)] 
					       + d45*(*field).u_comp[IND_2(n1f-6,k)]
					       + d46*(*field).u_comp[IND_2(n1f-7,k)]
					       )*h_1;
	    }
	    

/* 	    printf("at case 1"); */
}
	break;


	/* Take derivative of second coordinate */
	
	case 2:{ /* case 2: */
	    register int k;
	    for (k = n1i; k < n1f; k++){// for (k

		/* Value at 0 */

		(*dfield).u_comp[IND_2(k,n2i)] = (d00*(*field).u_comp[IND_2(k,n2i)] 
					    + d01*(*field).u_comp[IND_2(k,n2i+1)]
					    + d02*(*field).u_comp[IND_2(k,n2i+2)] 
					    + d03*(*field).u_comp[IND_2(k,n2i+3)]
					    )*h_2;
    
		/* Value at 1 */
		
		(*dfield).u_comp[IND_2(k,n2i+1)] = (d10*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d11*(*field).u_comp[IND_2(k,n2i+1)]
					      + d12*(*field).u_comp[IND_2(k,n2i+2)]
					      + d13*(*field).u_comp[IND_2(k,n2i+3)]
					      + d14*(*field).u_comp[IND_2(k,n2i+4)]
					      + d15*(*field).u_comp[IND_2(k,n2i+5)]
					      )*h_2;
		
		/* Value at 2 */
    
		(*dfield).u_comp[IND_2(k,n2i+2)] = (d20*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d21*(*field).u_comp[IND_2(k,n2i+1)]
					      + d22*(*field).u_comp[IND_2(k,n2i+2)]
					      + d23*(*field).u_comp[IND_2(k,n2i+3)] 
					      + d24*(*field).u_comp[IND_2(k,n2i+4)]
					      + d25*(*field).u_comp[IND_2(k,n2i+5)]
					      )*h_2;

		/* Value at 3 */

		(*dfield).u_comp[IND_2(k,n2i+3)] = (d30*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d31*(*field).u_comp[IND_2(k,n2i+1)]
					      + d32*(*field).u_comp[IND_2(k,n2i+2)]
					      + d33*(*field).u_comp[IND_2(k,n2i+3)]
					      + d34*(*field).u_comp[IND_2(k,n2i+4)] 
					      + d35*(*field).u_comp[IND_2(k,n2i+5)]
					      + d36*(*field).u_comp[IND_2(k,n2i+6)]
					      )*h_2;


		/* Value at 4 */

		(*dfield).u_comp[IND_2(k,n2i+4)] = (d40*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d41*(*field).u_comp[IND_2(k,n2i+1)]
					      + d42*(*field).u_comp[IND_2(k,n2i+2)]
					      + d43*(*field).u_comp[IND_2(k,n2i+3)]
					      + d44*(*field).u_comp[IND_2(k,n2i+4)] 
					      + d45*(*field).u_comp[IND_2(k,n2i+5)]
					      + d46*(*field).u_comp[IND_2(k,n2i+6)]
					      )*h_2;


		/* Intermediate Values 5---(n-6) */

		{register int j; 
		for (j = n2i+5; j < n2f-5; ++j){
		    (*dfield).u_comp[IND_2(k,j)] = (-d2*(*field).u_comp[IND_2(k,j-2)] 
					      - d1*(*field).u_comp[IND_2(k,j-1)]
					      + d1*(*field).u_comp[IND_2(k,j+1)] 
					      + d2*(*field).u_comp[IND_2(k,j+2)]
					      )*h_2;
		}
		}
		
		/* Value at n-1 */
		
		(*dfield).u_comp[IND_2(k,n2f-1)] = -(d00*(*field).u_comp[IND_2(k,n2f-1)] 
					       + d01*(*field).u_comp[IND_2(k,n2f-2)]
					       + d02*(*field).u_comp[IND_2(k,n2f-3)] 
					       + d03*(*field).u_comp[IND_2(k,n2f-4)]
					       )*h_2;
		
		/* Value at n-2 */
		
		(*dfield).u_comp[IND_2(k,n2f-2)] = -(d10*(*field).u_comp[IND_2(k,n2f-1)] 
					       + d11*(*field).u_comp[IND_2(k,n2f-2)]
					       + d12*(*field).u_comp[IND_2(k,n2f-3)]
					       + d13*(*field).u_comp[IND_2(k,n2f-4)]
					       + d14*(*field).u_comp[IND_2(k,n2f-5)]
					       + d15*(*field).u_comp[IND_2(k,n2f-6)]
					       )*h_2;
		
		/* Value at n-3 */
		
		(*dfield).u_comp[IND_2(k,n2f-3)] = -( d20*(*field).u_comp[IND_2(k,n2f-1)] 
						+ d21*(*field).u_comp[IND_2(k,n2f-2)]
						+ d22*(*field).u_comp[IND_2(k,n2f-3)]
						+ d23*(*field).u_comp[IND_2(k,n2f-4)]
						+ d24*(*field).u_comp[IND_2(k,n2f-5)]
						+ d25*(*field).u_comp[IND_2(k,n2f-6)]
						)*h_2;

		/* Value at n-4 */

		(*dfield).u_comp[IND_2(k,n2f-4)] = -(d30*(*field).u_comp[IND_2(k,n2f-1)]
					       + d31*(*field).u_comp[IND_2(k,n2f-2)]
					       + d32*(*field).u_comp[IND_2(k,n2f-3)]
					       + d33*(*field).u_comp[IND_2(k,n2f-4)]
					       + d34*(*field).u_comp[IND_2(k,n2f-5)] 
					       + d35*(*field).u_comp[IND_2(k,n2f-6)]
					       + d36*(*field).u_comp[IND_2(k,n2f-7)]
					       )*h_2;


		/* Value at n-5 */

		(*dfield).u_comp[IND_2(k,n2f-5)] = -(d40*(*field).u_comp[IND_2(k,n2f-1)]
					       + d41*(*field).u_comp[IND_2(k,n2f-2)]
					       + d42*(*field).u_comp[IND_2(k,n2f-3)]
					       + d43*(*field).u_comp[IND_2(k,n2f-4)]
					       + d44*(*field).u_comp[IND_2(k,n2f-5)] 
					       + d45*(*field).u_comp[IND_2(k,n2f-6)]
					       + d46*(*field).u_comp[IND_2(k,n2f-7)]
					       )*h_2;

	    }
	    

/* 	    printf("at case 2"); */
	}
	break;
	

	

	
	default:{
	    printf("wrong value given for coordinate\n");
	    exit(0);
	}
	
	printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
	printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
	}
    
}


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */





void deriv_strand_third_order_boundaries_sixth_interior_2d(struct GRID_PAR *grid,
							 struct field *field,
							 struct field *dfield,
							 int *coordinate)
{// void deriv_strand_third_order_boundaries_sixth_interior_2d
  
  int n1 = (*grid).n_gridpts_1;
  int n2 = (*grid).n_gridpts_2;
  
  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;

  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

#ifdef EXCISION
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");  
#endif

  int coord = (*coordinate);

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  /*----------------------------------------------------------------------------------------------- 
     I am writing down the closed form expressions (as fractions), for the coefficients of the 
     difference operator, the d's (called q's in Strand's paper), so that it is easier to debug, if ever needed. 
     The fractions are taken from Strand's paper, I checked the calculations and it seems to me that there 
     are no typos there, neither in the difference operator, nor in the scalar product.
    

     d00=-21600/13649, d01=81763/40947, d02=131/27298, d03=-9143/13649, d04=20539/81894

     d10=-81763/180195, d12=7357/36039, d13=30637/72078, d14=-2328/12013, d15=6611/360390

     d20=-131/54220, d21=-7357/16266, d23=645/2711, d24=11237/32532, d25=-3487/27110

     d30=9143/53590, d31=-30637/64308, d32=-645/5359, d34=-13733/32154, d35=-67/4660, d36=72/5359
     
     d40=-20539/236310, d41=2328/7877, d42=-11237/47262, d43=-13733/23631, d45=89387/118155,
     d46=-1296/7877, d47=144/7877

     d51=-6611/262806, d52=3487/43801, d53=1541/87602, d54=-89387/131403, d56=32400/43801, 
     d57=-6480/43801, d58=720/43801

  ----------------------------------------------------------------------------------------------*/
       //When computed with maple with 32 digits, these numbers are now

  
  FLOAT d00= -1.5825335189391164187852589933328;
  FLOAT d01= 1.9968007424231323418077026399980;
  FLOAT d02= 0.0047988863653014872884460400029306;
  FLOAT d03= -0.66986592424353432485896402666862;
  FLOAT d04= 0.25079981439421691454807434000049;
  
  FLOAT d10= -0.45374732928216654180193679069897;
  FLOAT d12= 0.20413995948833208468603457365632;
  FLOAT d13= 0.42505341435666916396126418602070;
  FLOAT d14= -0.19379006076750187297094813951552; 
  FLOAT d15= 0.018344016204667166125586170537473;

  FLOAT d20= -0.0024160826263371449649575802286979;
  FLOAT d21= -0.45229312676749047092093938276159;
  FLOAT d23= 0.23791958686831427517521209885651;
  FLOAT d24= 0.34541374646501905815812123447682;
  FLOAT d25= -0.12862412393950571744743637034305;

  FLOAT d30= 0.17061018846799776077626422840082;
  FLOAT d31= -0.47641039995023947253840890713442;
  FLOAT d32= -0.12035827579772345586863220750140;
  FLOAT d34= 0.42710082726876904895191889034024;
  FLOAT d35= -0.014377682403433476394849785407725;
  FLOAT d36= 0.013435342414629595073707781302482;

  FLOAT d40= -0.086915492361728238331005882104016;
  FLOAT d41= 0.29554398882823409927637425415767;
  FLOAT d42= -0.23775972239854428504929964876645;
  FLOAT d43= -0.58114341331302103169565401379544;
  FLOAT d45= 0.75652321103635055647243028225636;
  FLOAT d46= -0.16452964326520248825695061571664;
  FLOAT d47= 0.018281071473911387584105623968516;

  FLOAT d51= -0.025155437851495019139593464380570;
  FLOAT d52= 0.079610054564964270222141047008059;
  FLOAT d53= 0.017590922581676217437958037487729;
  FLOAT d54= -0.68025083141176381056749084876297;
  FLOAT d56= 0.73970913906075203762471176457159;
  FLOAT d57= -0.14794182781215040752494235291432;
  FLOAT d58= 0.016437980868016711947215816990480;

  FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
  FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);

  // these are for the intermediate values 
  FLOAT d1= 0.750000000000000000000000000000000;
  FLOAT d2= -0.15000000000000000000000000000000;
  FLOAT d3= 0.016666666666666666666666666666667; 

  switch(coord)
      { // switch(coord)

      /* Take derivative of first coordinate */

    case 1:{ /* Take derivative of first coordinate */
      {register int k;
      for (k = n2i; k < n2f; k++){
	  //---------------------->>>> AT AND CLOSE TO LEFT BOUNDARY 
	  /* Value at 0 */
	  (*dfield).u_comp[IND_2(n1i,k)] = (d00*(*field).u_comp[IND_2(n1i,k)] 
					      + d01*(*field).u_comp[IND_2(n1i+1,k)]
					      + d02*(*field).u_comp[IND_2(n1i+2,k)] 
					      + d03*(*field).u_comp[IND_2(n1i+3,k)]
					      + d04*(*field).u_comp[IND_2(n1i+4,k)]
					      )*h_1;
	  /* Value at 1 */
	  (*dfield).u_comp[IND_2(n1i+1,k)] = (d10*(*field).u_comp[IND_2(n1i+0,k)]
						+ d12*(*field).u_comp[IND_2(n1i+2,k)]
						+ d13*(*field).u_comp[IND_2(n1i+3,k)]
						+ d14*(*field).u_comp[IND_2(n1i+4,k)]
						+ d15*(*field).u_comp[IND_2(n1i+5,k)]
						)*h_1;
	  /* Value at 2 */
	  (*dfield).u_comp[IND_2(n1i+2,k)] = (d20*(*field).u_comp[IND_2(n1i+0,k)] 
						+ d21*(*field).u_comp[IND_2(n1i+1,k)]
						+ d23*(*field).u_comp[IND_2(n1i+3,k)] 
						+ d24*(*field).u_comp[IND_2(n1i+4,k)]
						+ d25*(*field).u_comp[IND_2(n1i+5,k)]
						)*h_1;
	  /* Value at 3 */
	  (*dfield).u_comp[IND_2(n1i+3,k)] = (d30*(*field).u_comp[IND_2(n1i+0,k)] 	
						+ d31*(*field).u_comp[IND_2(n1i+1,k)]
						+ d32*(*field).u_comp[IND_2(n1i+2,k)]
						+ d34*(*field).u_comp[IND_2(n1i+4,k)] 
						+ d35*(*field).u_comp[IND_2(n1i+5,k)]
						+ d36*(*field).u_comp[IND_2(n1i+6,k)]
						)*h_1;
	  /* Value at 4 */
	  (*dfield).u_comp[IND_2(n1i+4,k)] = (d40*(*field).u_comp[IND_2(n1i+0,k)] 	
						+ d41*(*field).u_comp[IND_2(n1i+1,k)]
						+ d42*(*field).u_comp[IND_2(n1i+2,k)]
						+ d43*(*field).u_comp[IND_2(n1i+3,k)] 
						+ d45*(*field).u_comp[IND_2(n1i+5,k)]
						+ d46*(*field).u_comp[IND_2(n1i+6,k)]
						+ d47*(*field).u_comp[IND_2(n1i+7,k)]
						)*h_1;
	  /* Value at 5 */
	  (*dfield).u_comp[IND_2(n1i+5,k)] = (d51*(*field).u_comp[IND_2(n1i+1,k)] 	
						+ d52*(*field).u_comp[IND_2(n1i+2,k)]
						+ d53*(*field).u_comp[IND_2(n1i+3,k)]
						+ d54*(*field).u_comp[IND_2(n1i+4,k)] 
						+ d56*(*field).u_comp[IND_2(n1i+6,k)]
						+ d57*(*field).u_comp[IND_2(n1i+7,k)]
						+ d58*(*field).u_comp[IND_2(n1i+8,k)]
						)*h_1;


	  //------------------>>  INTERMEDIATE VALUES 6-->(n-7)
	  {register int j; 
	  for (j = n1i+6; j < n1f-6; ++j){
	    (*dfield).u_comp[IND_2(j,k)] = (-d3*(*field).u_comp[IND_2(j-3,k)] 
					      -d2*(*field).u_comp[IND_2(j-2,k)] 
					      - d1*(*field).u_comp[IND_2(j-1,k)]
					      + d1*(*field).u_comp[IND_2(j+1,k)] 
					      + d2*(*field).u_comp[IND_2(j+2,k)]
					      + d3*(*field).u_comp[IND_2(j+3,k)]
					      )*h_1;
	  }
	  }


	  //-----------------------> AT AND CLOSE TO RIGHT BOUNDARY
	  /* Value at n-1 */
	  (*dfield).u_comp[IND_2(n1f-1,k)] = -(d00*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d01*(*field).u_comp[IND_2(n1f-2,k)]
					      + d02*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d03*(*field).u_comp[IND_2(n1f-4,k)]
					      + d04*(*field).u_comp[IND_2(n1f-5,k)]
					      )*h_1;
	  /* Value at n-2 */
	  (*dfield).u_comp[IND_2(n1f-2,k)] = -(d10*(*field).u_comp[IND_2(n1f-1,k)]
						+ d12*(*field).u_comp[IND_2(n1f-3,k)]
						+ d13*(*field).u_comp[IND_2(n1f-4,k)]
						+ d14*(*field).u_comp[IND_2(n1f-5,k)]
						+ d15*(*field).u_comp[IND_2(n1f-6,k)]
						)*h_1;
	  /* Value at n-3 */
	  (*dfield).u_comp[IND_2(n1f-3,k)] = -(d20*(*field).u_comp[IND_2(n1f-1,k)] 
						+ d21*(*field).u_comp[IND_2(n1f-2,k)]
						+ d23*(*field).u_comp[IND_2(n1f-4,k)] 
						+ d24*(*field).u_comp[IND_2(n1f-5,k)]
						+ d25*(*field).u_comp[IND_2(n1f-6,k)]
						)*h_1;
	  /* Value at n-4 */
	  (*dfield).u_comp[IND_2(n1f-4,k)] = -(d30*(*field).u_comp[IND_2(n1f-1,k)] 	
						+ d31*(*field).u_comp[IND_2(n1f-2,k)]
						+ d32*(*field).u_comp[IND_2(n1f-3,k)]
						+ d34*(*field).u_comp[IND_2(n1f-5,k)] 
						+ d35*(*field).u_comp[IND_2(n1f-6,k)]
						+ d36*(*field).u_comp[IND_2(n1f-7,k)]
						)*h_1;
	  /* Value at n-5 */
	  (*dfield).u_comp[IND_2(n1f-5,k)] = -(d40*(*field).u_comp[IND_2(n1f-1,k)] 	
						+ d41*(*field).u_comp[IND_2(n1f-2,k)]
						+ d42*(*field).u_comp[IND_2(n1f-3,k)]
						+ d43*(*field).u_comp[IND_2(n1f-4,k)] 
						+ d45*(*field).u_comp[IND_2(n1f-6,k)]
						+ d46*(*field).u_comp[IND_2(n1f-7,k)]
						+ d47*(*field).u_comp[IND_2(n1f-8,k)]
						)*h_1;
	  /* Value at n-6 */
	  (*dfield).u_comp[IND_2(n1f-6,k)] = -(d51*(*field).u_comp[IND_2(n1f-2,k)] 	
						+ d52*(*field).u_comp[IND_2(n1f-3,k)]
						+ d53*(*field).u_comp[IND_2(n1f-4,k)]
						+ d54*(*field).u_comp[IND_2(n1f-5,k)] 
						+ d56*(*field).u_comp[IND_2(n1f-7,k)]
						+ d57*(*field).u_comp[IND_2(n1f-8,k)]
						+ d58*(*field).u_comp[IND_2(n1f-9,k)]
						)*h_1;
	}
      }


#ifdef EXCISION
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
#ifdef DEBUG
      printf("at case 1");
#endif
    }
	    
	
      break;
	

      /* Take derivative of second coordinate */

    case 2:{ /* Take derivative of second coordinate */
	    
      {register int k;
      for (k = n1i; k < n1f; k++){// for (k = n1i; k < n1f; k++)

 //---------------------->>>> AT AND CLOSE TO LEFT BOUNDARY 
	  /* Value at 0 */
	  (*dfield).u_comp[IND_2(k,n2i)] = (d00*(*field).u_comp[IND_2(k,n2i)] 
					      + d01*(*field).u_comp[IND_2(k,n2i+1)]
					      + d02*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d03*(*field).u_comp[IND_2(k,n2i+3)]
					      + d04*(*field).u_comp[IND_2(k,n2i+4)]
					      )*h_2;
	  /* Value at 1 */
	  (*dfield).u_comp[IND_2(k,n2i+1)] = (d10*(*field).u_comp[IND_2(k,n2i+0)]
						+ d12*(*field).u_comp[IND_2(k,n2i+2)]
						+ d13*(*field).u_comp[IND_2(k,n2i+3)]
						+ d14*(*field).u_comp[IND_2(k,n2i+4)]
						+ d15*(*field).u_comp[IND_2(k,n2i+5)]
						)*h_2;
	  /* Value at 2 */
	  (*dfield).u_comp[IND_2(k,n2i+2)] = (d20*(*field).u_comp[IND_2(k,n2i+0)] 
						+ d21*(*field).u_comp[IND_2(k,n2i+1)]
						+ d23*(*field).u_comp[IND_2(k,n2i+3)] 
						+ d24*(*field).u_comp[IND_2(k,n2i+4)]
						+ d25*(*field).u_comp[IND_2(k,n2i+5)]
						)*h_2;
	  /* Value at 3 */
	  (*dfield).u_comp[IND_2(k,n2i+3)] = (d30*(*field).u_comp[IND_2(k,n2i+0)] 	
						+ d31*(*field).u_comp[IND_2(k,n2i+1)]
						+ d32*(*field).u_comp[IND_2(k,n2i+2)]
						+ d34*(*field).u_comp[IND_2(k,n2i+4)] 
						+ d35*(*field).u_comp[IND_2(k,n2i+5)]
						+ d36*(*field).u_comp[IND_2(k,n2i+6)]
						)*h_2;
	  /* Value at 4 */
	  (*dfield).u_comp[IND_2(k,n2i+4)] = (d40*(*field).u_comp[IND_2(k,n2i+0)] 	
						+ d41*(*field).u_comp[IND_2(k,n2i+1)]
						+ d42*(*field).u_comp[IND_2(k,n2i+2)]
						+ d43*(*field).u_comp[IND_2(k,n2i+3)] 
						+ d45*(*field).u_comp[IND_2(k,n2i+5)]
						+ d46*(*field).u_comp[IND_2(k,n2i+6)]
						+ d47*(*field).u_comp[IND_2(k,n2i+7)]
						)*h_2;
	  /* Value at 5 */
	  (*dfield).u_comp[IND_2(k,n2i+5)] = (d51*(*field).u_comp[IND_2(k,n2i+1)] 	
						+ d52*(*field).u_comp[IND_2(k,n2i+2)]
						+ d53*(*field).u_comp[IND_2(k,n2i+3)]
						+ d54*(*field).u_comp[IND_2(k,n2i+4)] 
						+ d56*(*field).u_comp[IND_2(k,n2i+6)]
						+ d57*(*field).u_comp[IND_2(k,n2i+7)]
						+ d58*(*field).u_comp[IND_2(k,n2i+8)]
						)*h_2;


	  //------------------>>  INTERMEDIATE VALUES 6-->(n-7)
	  {register int j; 
	  for (j = n2i+6; j < n2f-6; ++j){
	    (*dfield).u_comp[IND_2(k,j)] = (-d3*(*field).u_comp[IND_2(k,j-3)] 
					      -d2*(*field).u_comp[IND_2(k,j-2)] 
					      - d1*(*field).u_comp[IND_2(k,j-1)]
					      + d1*(*field).u_comp[IND_2(k,j+1)] 
					      + d2*(*field).u_comp[IND_2(k,j+2)]
					      + d3*(*field).u_comp[IND_2(k,j+3)]
					      )*h_2;
	  }}
	  


	  //-----------------------> AT AND CLOSE TO RIGHT BOUNDARY
	  /* Value at n-1 */
	  (*dfield).u_comp[IND_2(k,n2f-1)] = -(d00*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d01*(*field).u_comp[IND_2(k,n2f-2)]
					      + d02*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d03*(*field).u_comp[IND_2(k,n2f-4)]
					      + d04*(*field).u_comp[IND_2(k,n2f-5)]
					      )*h_2;
	  /* Value at n-2 */
	  (*dfield).u_comp[IND_2(k,n2f-2)] = -(d10*(*field).u_comp[IND_2(k,n2f-1)]
						+ d12*(*field).u_comp[IND_2(k,n2f-3)]
						+ d13*(*field).u_comp[IND_2(k,n2f-4)]
						+ d14*(*field).u_comp[IND_2(k,n2f-5)]
						+ d15*(*field).u_comp[IND_2(k,n2f-6)]
						)*h_2;
	  /* Value at n-3 */
	  (*dfield).u_comp[IND_2(k,n2f-3)] = -(d20*(*field).u_comp[IND_2(k,n2f-1)] 
						+ d21*(*field).u_comp[IND_2(k,n2f-2)]
						+ d23*(*field).u_comp[IND_2(k,n2f-4)] 
						+ d24*(*field).u_comp[IND_2(k,n2f-5)]
						+ d25*(*field).u_comp[IND_2(k,n2f-6)]
						)*h_2;
	  /* Value at n-4 */
	  (*dfield).u_comp[IND_2(k,n2f-4)] = -(d30*(*field).u_comp[IND_2(k,n2f-1)] 	
						+ d31*(*field).u_comp[IND_2(k,n2f-2)]
						+ d32*(*field).u_comp[IND_2(k,n2f-3)]
						+ d34*(*field).u_comp[IND_2(k,n2f-5)] 
						+ d35*(*field).u_comp[IND_2(k,n2f-6)]
						+ d36*(*field).u_comp[IND_2(k,n2f-7)]
						)*h_2;
	  /* Value at n-5 */
	  (*dfield).u_comp[IND_2(k,n2f-5)] = -(d40*(*field).u_comp[IND_2(k,n2f-1)] 	
						+ d41*(*field).u_comp[IND_2(k,n2f-2)]
						+ d42*(*field).u_comp[IND_2(k,n2f-3)]
						+ d43*(*field).u_comp[IND_2(k,n2f-4)] 
						+ d45*(*field).u_comp[IND_2(k,n2f-6)]
						+ d46*(*field).u_comp[IND_2(k,n2f-7)]
						+ d47*(*field).u_comp[IND_2(k,n2f-8)]
						)*h_2;
	  /* Value at n-6 */
	  (*dfield).u_comp[IND_2(k,n2f-6)] = -(d51*(*field).u_comp[IND_2(k,n2f-2)] 	
						+ d52*(*field).u_comp[IND_2(k,n2f-3)]
						+ d53*(*field).u_comp[IND_2(k,n2f-4)]
						+ d54*(*field).u_comp[IND_2(k,n2f-5)] 
						+ d56*(*field).u_comp[IND_2(k,n2f-7)]
						+ d57*(*field).u_comp[IND_2(k,n2f-8)]
						+ d58*(*field).u_comp[IND_2(k,n2f-9)]
						)*h_2;

	  
	}
      }
      


#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
#ifdef DEBUG
      printf("at case 2");
#endif
    }

	
      break;

	
    default:{
      printf("wrong value given for coordinate\n");
      exit(0);
    }

      }
}



/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */







/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


//-----------------------------------------------------------------------------------------

void deriv_strand_fourth_order_boundaries_eight_interior_2d(struct GRID_PAR *grid,
							 struct field *field,
							 struct field *dfield,
							 int *coordinate)
{// deriv_strand_fourth_order_boundaries_eight_interior_2d
  

  
  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;

  if (n_gridpts_1<16||n_gridpts_2<16){
    printf("***************************************************************************************\n");
    printf("error: you are not using enough points for SBP to hold, you need at least 16 in each direction \n");
    printf("and you are using --------> %d %d  \n", n_gridpts_1, n_gridpts_2);
    printf("***************************************************************************************\n");
    exit(1); }
  
  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

#ifdef EXCISION
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");  
#endif

  int coord = (*coordinate);

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  /*----------------------------------------------------------------------------------------------- 
    The coeffs from Strand's paper for this case are Ok, I checked the calculations and it seems to me that there 
     are no typos there, neither in the difference operator, nor in the scalar product.
     ----------------------------------------------------------------------------------------------*/
  //When computed with maple with 32 digits, these numbers are now

#ifdef NON_OPTIMIZED_EIGHT_ORDER_OPERATOR
  FLOAT d00=-1.6955436044318985087498556542484;
  FLOAT d01=2.0610513554928258770826116045752;
  FLOAT d02=0.87789728901434824583477679084963;
  FLOAT d03=-2.5445639556810149125014434575163;
  FLOAT d04=1.6889486445071741229173883954248;
  FLOAT d05=-0.38778972890143482458347767908496;
  FLOAT d06= 0.0;
  FLOAT d07= 0.0;
  FLOAT d08= 0.0;
  FLOAT d09= 0.0;
  FLOAT d010= 0.0;
  FLOAT d011= 0.0;

  FLOAT d10=-0.39835918734972926809517745906661;
  FLOAT d11= 0.0;
  FLOAT d12=-0.44127885642072764523900478066757;
  FLOAT d13=1.8989658393441626095262349706691;
  FLOAT d14=-1.5738365692621829357170143420027;
  FLOAT d15=0.60604617027316423238240764906811;
  FLOAT d16=-0.091537396584686992857446038000302;
  FLOAT d17= 0.0;
  FLOAT d18= 0.0;
  FLOAT d19= 0.0;
  FLOAT d110= 0.0;
  FLOAT d111= 0.0;

  FLOAT d20=-1.0055577090642838920860210778808;
  FLOAT d21=2.6151125596961539353524325346492;
  FLOAT d22= 0.0;
  FLOAT d23=-10.448835244859039473942416960576;
  FLOAT d24=16.854276269429682078398421541815;
  FLOAT d25=-10.479793844227156688020808246231;
  FLOAT d26=2.5403284648744886691035720465902;
  FLOAT d27=-0.075530495849844628805179838367064;
  FLOAT d28= 0.0;
  FLOAT d29= 0.0;
  FLOAT d210= 0.0;
  FLOAT d211= 0.0;

  FLOAT d30=0.41730852995727528198434306408336;
  FLOAT d31=-1.6112948490126556929673253900322;
  FLOAT d32=1.4960581706703734383449080548465;
  FLOAT d33= 0.0;
  FLOAT d34=-1.6738166082885887268879133240715;
  FLOAT d35=1.9668117242490862700533359284822;
  FLOAT d36=-0.64585509260926636725392126737262;
  FLOAT d37=0.050788125033775796726572934064281;
  FLOAT d38= 0.0;
  FLOAT d39= 0.0;
  FLOAT d310= 0.0;
  FLOAT d311= 0.0;

  FLOAT d40=-1.2067978767157806109219504339143;
  FLOAT d41=5.8182409265274916785465083281307;
  FLOAT d42=-10.513925845304986417473772536922;
  FLOAT d43=7.2925946575767793888363987219850;
  FLOAT d44= 0.0;
  FLOAT d45=-3.0385143353798934097204370001146;
  FLOAT d46=2.0057652009112144592418490931814;
  FLOAT d47=-0.34870908370292399111170040001164;
  FLOAT d48=-0.0086536439119010973968957723343805;
  FLOAT d49= 0.0;
  FLOAT d410= 0.0;
  FLOAT d411= 0.0;

  FLOAT d50=0.089446187550400088107295882497877;
  FLOAT d51=-0.72324463225750363106445242685870;
  FLOAT d52=2.1103523865873885758072494834961;
  FLOAT d53=-2.7662054325681783318765966175687;
  FLOAT d54=0.98086385026111691291085453864196;
  FLOAT d55= 0.0;
  FLOAT d56=0.30802293704705996043280348352271;
  FLOAT d57=-0.026238992169135626268229544755535;
  FLOAT d58=0.029797181295285022842565739061272;
  FLOAT d59=-0.0027934857464329708914905380369943;
  FLOAT d510= 0.0;
  FLOAT d511= 0.0;

  FLOAT d60= 0.0;
  FLOAT d61=0.15126303740835199995224529205735;
  FLOAT d62=-0.70834831886017471258878769201994;
  FLOAT d63=1.2577996869081960509841200672744;
  FLOAT d64=-0.89656603854302375846464634405377;
  FLOAT d65=-0.42651843804299217071512177699248;
  FLOAT d66= 0.0;
  FLOAT d67=0.80159349003841766765011518336067;
  FLOAT d68=-0.21661535522787203529072916962025;
  FLOAT d69=0.041260067662451816245853175165761;
  FLOAT d610=-0.0038681313433548577730487351717901;
  FLOAT d611= 0.0;

  FLOAT d70= 0.0;
  FLOAT d71= 0.0;
  FLOAT d72=0.019265719907611002302054245220619;
  FLOAT d73=-0.090478311526126167237971095382715;
  FLOAT d74=0.14258418768974005892265577479665;
  FLOAT d75=0.033235928479719164073418453370319;
  FLOAT d76=-0.73326354624003548802568409455577;
  FLOAT d77= 0.0;
  FLOAT d78=0.79260196355547737511601116983528;
  FLOAT d79=-0.19815049088886934377900279245882;
  FLOAT d710=0.037742950645498922624571960468347;
  FLOAT d711=-0.0035384016230155239960536212939075;
#endif

  // these are the coeffs for the improved operator that we choose not by minimizing the bandwith, but 
  // by minimizing the maximum imaginary eigenvalue.

#ifdef OPTIMIZED_EIGHT_ORDER_OPERATOR
  FLOAT d00= -1.6955436044318985087498556542484;
  FLOAT d01= 2.252142137678813514633822362278;
  FLOAT d02= -0.05888015064022764242837280119;
  FLOAT d03= -0.72703849019795003445385686285;
  FLOAT d04= -0.03519446459907925766567721687;
  FLOAT d05= 0.3808994692748803682435341446956;
  FLOAT d06= -0.097708425809176140086689775338;
  FLOAT d07= -0.018676471275362299492904196473;
  FLOAT d08= 0.;
  FLOAT d09= 0.;
  FLOAT d010=0.;
  FLOAT d011= 0.;

  FLOAT d10= -0.4352931378302752275823881078793;
  FLOAT d11= 0.;
  FLOAT d12= 0.099854313212144418846805809049;
  FLOAT d13= 0.48598825799891087114955473054;
  FLOAT d14= -0.04056967339078804101602657017;
  FLOAT d15= -0.1516077655067655130871564977011;
  FLOAT d16= 0.028751917941456163466861493696;
  FLOAT d17= 0.012876087575317328222349142474;
  FLOAT d18= 0.;
  FLOAT d19= 0.;
  FLOAT d110= 0.;
  FLOAT d111= 0.;

  FLOAT d20= 0.06744227386055814358882611857;
  FLOAT d21= -0.59175794358011014241446301559;
  FLOAT d22= 0.;
  FLOAT d23= 0.09922290191545044394136441726;
  FLOAT d24= 1.24445434548753271358155515947;
  FLOAT d25= -1.368619045325369683541255238187;
  FLOAT d26= 0.68543519098309329993738259237;
  FLOAT d27= -0.13617772334115477509341003392;
  FLOAT d28= 0.;
  FLOAT d29= 0.;
  FLOAT d210= 0.;
  FLOAT d211= 0.;

  FLOAT d30= 0.119234324171531019965386581760;
  FLOAT d31= -0.41236675277145806241997120304;
  FLOAT d32= -0.01420667755300961916445912985;
  FLOAT d33= 0.;
  FLOAT d34= -0.117132848377661592077848335517;
  FLOAT d35= 0.6750458046328004849045861337825;
  FLOAT d36= -0.288099428939536862991210465807;
  FLOAT d37= 0.037525578837334631783516418677;
  FLOAT d38= 0.;
  FLOAT d39= 0.;
  FLOAT d310= 0.;
  FLOAT d311= 0.;

  FLOAT d40= 0.025147363295176347088194768037;
  FLOAT d41= 0.14998007970344354042651691946;
  FLOAT d42= -0.77630747812384192410033152270;
  FLOAT d43= 0.51033212364828546341398271274;
  FLOAT d44= 0.;
  FLOAT d45= -0.568060141489748837333529197702;
  FLOAT d46= 0.891763513806768671939424492616;
  FLOAT d47= -0.224201816928182164037362400114;
  FLOAT d48= -0.0086536439119010973968957723343805;
  FLOAT d49= 0.;
  FLOAT d410= 0.;
  FLOAT d411= 0.;

  FLOAT d50= -0.0878569049859194092294958469532;
  FLOAT d51= 0.1809259887937250033140311227933;
  FLOAT d52= 0.275603557814485387547649301144;
  FLOAT d53= -0.949412365700909499468138419805;
  FLOAT d54= 0.1833756882676831850910631088522;
  FLOAT d55= 0.;
  FLOAT d56= 0.42315722086966643064298670184389;
  FLOAT d57= -0.052796880607583149849171168899192;
  FLOAT d58= 0.029797181295285022842565739061272;
  FLOAT d59= -0.0027934857464329708914905380369943;
  FLOAT d510= 0.;
  FLOAT d511= 0.;

  FLOAT d60= 0.031207020141045519688121961545;
  FLOAT d61= -0.047511755865994320174424070761;
  FLOAT d62= -0.191127593117949651614482212113;
  FLOAT d63= 0.561072252375672358928364455971;
  FLOAT d64= -0.398613397284846386535560095569;
  FLOAT d65= -0.58594453589139385546142240382277;
  FLOAT d66= 0.;
  FLOAT d67= 0.81014142855224141198732709437972;
  FLOAT d68= -0.21661535522787203529072916962025;
  FLOAT d69= 0.041260067662451816245853175165761;
  FLOAT d610= -0.0038681313433548577730487351717901;
  FLOAT d611= 0.;

  FLOAT d70= 0.0054565862264050490869367571164;
  FLOAT d71= -0.019463641447689387726897436343;
  FLOAT d72= 0.034735133749982204632490070184;
  FLOAT d73= -0.066851276946818081029475174146;
  FLOAT d74= 0.091674221978406727279476067978;
  FLOAT d75= 0.066875790674993403525413442454852;
  FLOAT d76= -0.74108283592437134573347044379599;
  FLOAT d77= 0.;
  FLOAT d78= 0.79260196355547737511601116983528;
  FLOAT d79= -0.19815049088886934377900279245882;
  FLOAT d710= 0.037742950645498922624571960468347;
  FLOAT d711= -0.0035384016230155239960536212939075;
#endif



  FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
  FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);
  // these are for the intermediate values 
  FLOAT d1= 0.800000000000000000000000000000000;
  FLOAT d2= -0.2000000000000000000000000000000;
  FLOAT d3= 0.038095238095238095238095238095238;
  FLOAT d4= -0.0035714285714285714285714285714286;


  switch(coord)
      {// switch(coord)

      /* Take derivative of first coordinate */

    case 1:{// case 1:
      {register int k;
      for (k = n2i; k < n2f; k++){// for k
	  //---------------------->>>> AT AND CLOSE TO LEFT BOUNDARY 
	  /* Value at 0 */
	  (*dfield).u_comp[IND_2(n1i,k)] = (d00*(*field).u_comp[IND_2(n1i+0,k)] 
					      + d01*(*field).u_comp[IND_2(n1i+1,k)]
					      + d02*(*field).u_comp[IND_2(n1i+2,k)] 
					      + d03*(*field).u_comp[IND_2(n1i+3,k)]
					      + d04*(*field).u_comp[IND_2(n1i+4,k)]
					      + d05*(*field).u_comp[IND_2(n1i+5,k)]
					      + d06*(*field).u_comp[IND_2(n1i+6,k)] 
					      + d07*(*field).u_comp[IND_2(n1i+7,k)]
					      + d08*(*field).u_comp[IND_2(n1i+8,k)] 
					      + d09*(*field).u_comp[IND_2(n1i+9,k)]
					      + d010*(*field).u_comp[IND_2(n1i+10,k)]
					      + d011*(*field).u_comp[IND_2(n1i+11,k)]
					      )*h_1;
	  /* Value at 1 */
	  (*dfield).u_comp[IND_2(n1i+1,k)] = (d10*(*field).u_comp[IND_2(n1i+0,k)]
						+ d11*(*field).u_comp[IND_2(n1i+1,k)]
						+ d12*(*field).u_comp[IND_2(n1i+2,k)]
						+ d13*(*field).u_comp[IND_2(n1i+3,k)]
						+ d14*(*field).u_comp[IND_2(n1i+4,k)]
						+ d15*(*field).u_comp[IND_2(n1i+5,k)]
						+ d16*(*field).u_comp[IND_2(n1i+6,k)]
						+ d17*(*field).u_comp[IND_2(n1i+7,k)]
						+ d18*(*field).u_comp[IND_2(n1i+8,k)] 
						+ d19*(*field).u_comp[IND_2(n1i+9,k)]
						+ d110*(*field).u_comp[IND_2(n1i+10,k)]
						+ d111*(*field).u_comp[IND_2(n1i+11,k)]
						)*h_1;
	  /* Value at 2 */
	  (*dfield).u_comp[IND_2(n1i+2,k)] = (d20*(*field).u_comp[IND_2(n1i+0,k)] 
						+ d21*(*field).u_comp[IND_2(n1i+1,k)]
						+ d22*(*field).u_comp[IND_2(n1i+2,k)]
						+ d23*(*field).u_comp[IND_2(n1i+3,k)] 
						+ d24*(*field).u_comp[IND_2(n1i+4,k)]
						+ d25*(*field).u_comp[IND_2(n1i+5,k)]
						+ d26*(*field).u_comp[IND_2(n1i+6,k)]
						+ d27*(*field).u_comp[IND_2(n1i+7,k)]
						+ d28*(*field).u_comp[IND_2(n1i+8,k)] 
						+ d29*(*field).u_comp[IND_2(n1i+9,k)]
						+ d210*(*field).u_comp[IND_2(n1i+10,k)]
						+ d211*(*field).u_comp[IND_2(n1i+11,k)]
						)*h_1;
	  /* Value at 3 */
	  (*dfield).u_comp[IND_2(n1i+3,k)] = (d30*(*field).u_comp[IND_2(n1i+0,k)] 	
						+ d31*(*field).u_comp[IND_2(n1i+1,k)]
						+ d32*(*field).u_comp[IND_2(n1i+2,k)]
						+ d33*(*field).u_comp[IND_2(n1i+3,k)]
						+ d34*(*field).u_comp[IND_2(n1i+4,k)] 
						+ d35*(*field).u_comp[IND_2(n1i+5,k)]
						+ d36*(*field).u_comp[IND_2(n1i+6,k)]
						+ d37*(*field).u_comp[IND_2(n1i+7,k)]
						+ d38*(*field).u_comp[IND_2(n1i+8,k)] 
						+ d39*(*field).u_comp[IND_2(n1i+9,k)]
						+ d310*(*field).u_comp[IND_2(n1i+10,k)]
						+ d311*(*field).u_comp[IND_2(n1i+11,k)]
						)*h_1;
	  /* Value at 4 */
	  (*dfield).u_comp[IND_2(n1i+4,k)] = (d40*(*field).u_comp[IND_2(n1i+0,k)] 	
						+ d41*(*field).u_comp[IND_2(n1i+1,k)]
						+ d42*(*field).u_comp[IND_2(n1i+2,k)]
						+ d43*(*field).u_comp[IND_2(n1i+3,k)]
						+ d44*(*field).u_comp[IND_2(n1i+4,k)]
						+ d45*(*field).u_comp[IND_2(n1i+5,k)]
						+ d46*(*field).u_comp[IND_2(n1i+6,k)]
						+ d47*(*field).u_comp[IND_2(n1i+7,k)]
						+ d48*(*field).u_comp[IND_2(n1i+8,k)]
						+ d49*(*field).u_comp[IND_2(n1i+9,k)]
						+ d410*(*field).u_comp[IND_2(n1i+10,k)]
						+ d411*(*field).u_comp[IND_2(n1i+11,k)]
						)*h_1;
	  /* Value at 5 */
	  (*dfield).u_comp[IND_2(n1i+5,k)] = (d50*(*field).u_comp[IND_2(n1i+0,k)] 
						+ d51*(*field).u_comp[IND_2(n1i+1,k)] 	
						+ d52*(*field).u_comp[IND_2(n1i+2,k)]
						+ d53*(*field).u_comp[IND_2(n1i+3,k)]
						+ d54*(*field).u_comp[IND_2(n1i+4,k)] 
						+ d55*(*field).u_comp[IND_2(n1i+5,k)] 
						+ d56*(*field).u_comp[IND_2(n1i+6,k)]
						+ d57*(*field).u_comp[IND_2(n1i+7,k)]
						+ d58*(*field).u_comp[IND_2(n1i+8,k)]
						+ d59*(*field).u_comp[IND_2(n1i+9,k)]
						+ d510*(*field).u_comp[IND_2(n1i+10,k)]
						+ d511*(*field).u_comp[IND_2(n1i+11,k)]
						)*h_1;

	  /* Value at 6 */
	  (*dfield).u_comp[IND_2(n1i+6,k)] = (d60*(*field).u_comp[IND_2(n1i+0,k)] 
						+ d61*(*field).u_comp[IND_2(n1i+1,k)] 	
						+ d62*(*field).u_comp[IND_2(n1i+2,k)]
						+ d63*(*field).u_comp[IND_2(n1i+3,k)]
						+ d64*(*field).u_comp[IND_2(n1i+4,k)] 
						+ d65*(*field).u_comp[IND_2(n1i+5,k)] 
						+ d66*(*field).u_comp[IND_2(n1i+6,k)] 
						+ d67*(*field).u_comp[IND_2(n1i+7,k)]
						+ d68*(*field).u_comp[IND_2(n1i+8,k)]
						+ d69*(*field).u_comp[IND_2(n1i+9,k)]
						+ d610*(*field).u_comp[IND_2(n1i+10,k)]
						+ d611*(*field).u_comp[IND_2(n1i+11,k)] 
						)*h_1;

	  /* Value at 7 */
	  (*dfield).u_comp[IND_2(n1i+7,k)] = (d70*(*field).u_comp[IND_2(n1i+0,k)] 
						+ d71*(*field).u_comp[IND_2(n1i+1,k)] 
						+ d72*(*field).u_comp[IND_2(n1i+2,k)]
						+ d73*(*field).u_comp[IND_2(n1i+3,k)]
						+ d74*(*field).u_comp[IND_2(n1i+4,k)] 
						+ d75*(*field).u_comp[IND_2(n1i+5,k)] 
						+ d76*(*field).u_comp[IND_2(n1i+6,k)] 
						+ d77*(*field).u_comp[IND_2(n1i+7,k)] 
						+ d78*(*field).u_comp[IND_2(n1i+8,k)]
						+ d79*(*field).u_comp[IND_2(n1i+9,k)]
						+ d710*(*field).u_comp[IND_2(n1i+10,k)]
						+ d711*(*field).u_comp[IND_2(n1i+11,k)]
						)*h_1;


	  //------------------>>  INTERMEDIATE VALUES 6-->(n-7)
	  {register int j; 
	  for (j = n1i+8; j < n1f-8; ++j){
	    (*dfield).u_comp[IND_2(j,k)] = (-d4*(*field).u_comp[IND_2(j-4,k)] 
					      -d3*(*field).u_comp[IND_2(j-3,k)] 
					      -d2*(*field).u_comp[IND_2(j-2,k)] 
					      - d1*(*field).u_comp[IND_2(j-1,k)]
					      + d1*(*field).u_comp[IND_2(j+1,k)] 
					      + d2*(*field).u_comp[IND_2(j+2,k)]
					      + d3*(*field).u_comp[IND_2(j+3,k)]
					      + d4*(*field).u_comp[IND_2(j+4,k)]
					      )*h_1;
	  }
	  }


	  //-----------------------> AT AND CLOSE TO RIGHT BOUNDARY
	  
	  /* Value at n-1 */
	  (*dfield).u_comp[IND_2(n1f-1,k)] = -(d00*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d01*(*field).u_comp[IND_2(n1f-2,k)]
					      + d02*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d03*(*field).u_comp[IND_2(n1f-4,k)]
					      + d04*(*field).u_comp[IND_2(n1f-5,k)]
					      + d05*(*field).u_comp[IND_2(n1f-6,k)]
					      + d06*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d07*(*field).u_comp[IND_2(n1f-8,k)]
					      + d08*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d09*(*field).u_comp[IND_2(n1f-10,k)]
					      + d010*(*field).u_comp[IND_2(n1f-11,k)]
					      + d011*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;
	  /* Value at n-2 */
	  (*dfield).u_comp[IND_2(n1f-2,k)] = -(d10*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d11*(*field).u_comp[IND_2(n1f-2,k)]
					      + d12*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d13*(*field).u_comp[IND_2(n1f-4,k)]
					      + d14*(*field).u_comp[IND_2(n1f-5,k)]
					      + d15*(*field).u_comp[IND_2(n1f-6,k)]
					      + d16*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d17*(*field).u_comp[IND_2(n1f-8,k)]
					      + d18*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d19*(*field).u_comp[IND_2(n1f-10,k)]
					      + d110*(*field).u_comp[IND_2(n1f-11,k)]
					      + d111*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;

	  /* Value at n-3 */
	  (*dfield).u_comp[IND_2(n1f-3,k)] = -(d20*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d21*(*field).u_comp[IND_2(n1f-2,k)]
					      + d22*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d23*(*field).u_comp[IND_2(n1f-4,k)]
					      + d24*(*field).u_comp[IND_2(n1f-5,k)]
					      + d25*(*field).u_comp[IND_2(n1f-6,k)]
					      + d26*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d27*(*field).u_comp[IND_2(n1f-8,k)]
					      + d28*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d29*(*field).u_comp[IND_2(n1f-10,k)]
					      + d210*(*field).u_comp[IND_2(n1f-11,k)]
					      + d211*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;

	  /* Value at n-4 */
	  (*dfield).u_comp[IND_2(n1f-4,k)] = -(d30*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d31*(*field).u_comp[IND_2(n1f-2,k)]
					      + d32*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d33*(*field).u_comp[IND_2(n1f-4,k)]
					      + d34*(*field).u_comp[IND_2(n1f-5,k)]
					      + d35*(*field).u_comp[IND_2(n1f-6,k)]
					      + d36*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d37*(*field).u_comp[IND_2(n1f-8,k)]
					      + d38*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d39*(*field).u_comp[IND_2(n1f-10,k)]
					      + d310*(*field).u_comp[IND_2(n1f-11,k)]
					      + d311*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;

	    /* Value at n-5 */
	    (*dfield).u_comp[IND_2(n1f-5,k)] = -(d40*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d41*(*field).u_comp[IND_2(n1f-2,k)]
					      + d42*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d43*(*field).u_comp[IND_2(n1f-4,k)]
					      + d44*(*field).u_comp[IND_2(n1f-5,k)]
					      + d45*(*field).u_comp[IND_2(n1f-6,k)]
					      + d46*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d47*(*field).u_comp[IND_2(n1f-8,k)]
					      + d48*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d49*(*field).u_comp[IND_2(n1f-10,k)]
					      + d410*(*field).u_comp[IND_2(n1f-11,k)]
					      + d411*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;
	  /* Value at n-6 */
	  (*dfield).u_comp[IND_2(n1f-6,k)] = -(d50*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d51*(*field).u_comp[IND_2(n1f-2,k)]
					      + d52*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d53*(*field).u_comp[IND_2(n1f-4,k)]
					      + d54*(*field).u_comp[IND_2(n1f-5,k)]
					      + d55*(*field).u_comp[IND_2(n1f-6,k)]
					      + d56*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d57*(*field).u_comp[IND_2(n1f-8,k)]
					      + d58*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d59*(*field).u_comp[IND_2(n1f-10,k)]
					      + d510*(*field).u_comp[IND_2(n1f-11,k)]
					      + d511*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;
	  /* Value at n-7 */
	  (*dfield).u_comp[IND_2(n1f-7,k)] = -(d60*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d61*(*field).u_comp[IND_2(n1f-2,k)]
					      + d62*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d63*(*field).u_comp[IND_2(n1f-4,k)]
					      + d64*(*field).u_comp[IND_2(n1f-5,k)]
					      + d65*(*field).u_comp[IND_2(n1f-6,k)]
					      + d66*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d67*(*field).u_comp[IND_2(n1f-8,k)]
					      + d68*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d69*(*field).u_comp[IND_2(n1f-10,k)]
					      + d610*(*field).u_comp[IND_2(n1f-11,k)]
					      + d611*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;

	  /* Value at n-8 */
	  (*dfield).u_comp[IND_2(n1f-8,k)] =  -(d70*(*field).u_comp[IND_2(n1f-1,k)] 
					      + d71*(*field).u_comp[IND_2(n1f-2,k)]
					      + d72*(*field).u_comp[IND_2(n1f-3,k)] 
					      + d73*(*field).u_comp[IND_2(n1f-4,k)]
					      + d74*(*field).u_comp[IND_2(n1f-5,k)]
					      + d75*(*field).u_comp[IND_2(n1f-6,k)]
					      + d76*(*field).u_comp[IND_2(n1f-7,k)] 
					      + d77*(*field).u_comp[IND_2(n1f-8,k)]
					      + d78*(*field).u_comp[IND_2(n1f-9,k)] 
					      + d79*(*field).u_comp[IND_2(n1f-10,k)]
					      + d710*(*field).u_comp[IND_2(n1f-11,k)]
					      + d711*(*field).u_comp[IND_2(n1f-12,k)]
					      )*h_1;


	
      }
      }


#ifdef EXCISION
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
#ifdef DEBUG
      printf("at case 1");
#endif
    }
	    
	
      break;
	

      /* Take derivative of second coordinate */

    case 2:{// case 2:
	    
      {register int k;
      for (k = n1i; k < n1f; k++){// for (k =
	 
	  //---------------------->>>> AT AND CLOSE TO LEFT BOUNDARY 
	  /* Value at 0 */
	  (*dfield).u_comp[IND_2(k,n2i)] = (d00*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d01*(*field).u_comp[IND_2(k,n2i+1)]
					      + d02*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d03*(*field).u_comp[IND_2(k,n2i+3)]
					      + d04*(*field).u_comp[IND_2(k,n2i+4)]
					      + d05*(*field).u_comp[IND_2(k,n2i+5)]
					      + d06*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d07*(*field).u_comp[IND_2(k,n2i+7)]
					      + d08*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d09*(*field).u_comp[IND_2(k,n2i+9)]
					      + d010*(*field).u_comp[IND_2(k,n2i+10)]
					      + d011*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;
	  /* Value at 1 */
	  (*dfield).u_comp[IND_2(k,n2i+1)] = (d10*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d11*(*field).u_comp[IND_2(k,n2i+1)]
					      + d12*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d13*(*field).u_comp[IND_2(k,n2i+3)]
					      + d14*(*field).u_comp[IND_2(k,n2i+4)]
					      + d15*(*field).u_comp[IND_2(k,n2i+5)]
					      + d16*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d17*(*field).u_comp[IND_2(k,n2i+7)]
					      + d18*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d19*(*field).u_comp[IND_2(k,n2i+9)]
					      + d110*(*field).u_comp[IND_2(k,n2i+10)]
					      + d111*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;
	  /* Value at 2 */
	  (*dfield).u_comp[IND_2(k,n2i+2)] = (d20*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d21*(*field).u_comp[IND_2(k,n2i+1)]
					      + d22*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d23*(*field).u_comp[IND_2(k,n2i+3)]
					      + d24*(*field).u_comp[IND_2(k,n2i+4)]
					      + d25*(*field).u_comp[IND_2(k,n2i+5)]
					      + d26*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d27*(*field).u_comp[IND_2(k,n2i+7)]
					      + d28*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d29*(*field).u_comp[IND_2(k,n2i+9)]
					      + d210*(*field).u_comp[IND_2(k,n2i+10)]
					      + d211*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;
	  /* Value at 3 */
	  (*dfield).u_comp[IND_2(k,n2i+3)] = (d30*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d31*(*field).u_comp[IND_2(k,n2i+1)]
					      + d32*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d33*(*field).u_comp[IND_2(k,n2i+3)]
					      + d34*(*field).u_comp[IND_2(k,n2i+4)]
					      + d35*(*field).u_comp[IND_2(k,n2i+5)]
					      + d36*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d37*(*field).u_comp[IND_2(k,n2i+7)]
					      + d38*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d39*(*field).u_comp[IND_2(k,n2i+9)]
					      + d310*(*field).u_comp[IND_2(k,n2i+10)]
					      + d311*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;
	  /* Value at 4 */
	  (*dfield).u_comp[IND_2(k,n2i+4)] = (d40*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d41*(*field).u_comp[IND_2(k,n2i+1)]
					      + d42*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d43*(*field).u_comp[IND_2(k,n2i+3)]
					      + d44*(*field).u_comp[IND_2(k,n2i+4)]
					      + d45*(*field).u_comp[IND_2(k,n2i+5)]
					      + d46*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d47*(*field).u_comp[IND_2(k,n2i+7)]
					      + d48*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d49*(*field).u_comp[IND_2(k,n2i+9)]
					      + d410*(*field).u_comp[IND_2(k,n2i+10)]
					      + d411*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;
	  /* Value at 5 */
	  (*dfield).u_comp[IND_2(k,n2i+5)] = (d50*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d51*(*field).u_comp[IND_2(k,n2i+1)]
					      + d52*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d53*(*field).u_comp[IND_2(k,n2i+3)]
					      + d54*(*field).u_comp[IND_2(k,n2i+4)]
					      + d55*(*field).u_comp[IND_2(k,n2i+5)]
					      + d56*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d57*(*field).u_comp[IND_2(k,n2i+7)]
					      + d58*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d59*(*field).u_comp[IND_2(k,n2i+9)]
					      + d510*(*field).u_comp[IND_2(k,n2i+10)]
					      + d511*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;

	  /* Value at 6 */
	  (*dfield).u_comp[IND_2(k,n2i+6)] = (d60*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d61*(*field).u_comp[IND_2(k,n2i+1)]
					      + d62*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d63*(*field).u_comp[IND_2(k,n2i+3)]
					      + d64*(*field).u_comp[IND_2(k,n2i+4)]
					      + d65*(*field).u_comp[IND_2(k,n2i+5)]
					      + d66*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d67*(*field).u_comp[IND_2(k,n2i+7)]
					      + d68*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d69*(*field).u_comp[IND_2(k,n2i+9)]
					      + d610*(*field).u_comp[IND_2(k,n2i+10)]
					      + d611*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;

	  /* Value at 7 */
	  (*dfield).u_comp[IND_2(k,n2i+7)] = (d70*(*field).u_comp[IND_2(k,n2i+0)] 
					      + d71*(*field).u_comp[IND_2(k,n2i+1)]
					      + d72*(*field).u_comp[IND_2(k,n2i+2)] 
					      + d73*(*field).u_comp[IND_2(k,n2i+3)]
					      + d74*(*field).u_comp[IND_2(k,n2i+4)]
					      + d75*(*field).u_comp[IND_2(k,n2i+5)]
					      + d76*(*field).u_comp[IND_2(k,n2i+6)] 
					      + d77*(*field).u_comp[IND_2(k,n2i+7)]
					      + d78*(*field).u_comp[IND_2(k,n2i+8)] 
					      + d79*(*field).u_comp[IND_2(k,n2i+9)]
					      + d710*(*field).u_comp[IND_2(k,n2i+10)]
					      + d711*(*field).u_comp[IND_2(k,n2i+11)]
					      )*h_2;



	  //------------------>>  INTERMEDIATE VALUES 6-->(n-7)
	  {register int j; 
	  for (j = n2i+8; j < n2f-8; ++j){
	    (*dfield).u_comp[IND_2(k,j)] = (-d4*(*field).u_comp[IND_2(k,j-4)] 
					      -d3*(*field).u_comp[IND_2(k,j-3)] 
					      -d2*(*field).u_comp[IND_2(k,j-2)] 
					      - d1*(*field).u_comp[IND_2(k,j-1)]
					      + d1*(*field).u_comp[IND_2(k,j+1)] 
					      + d2*(*field).u_comp[IND_2(k,j+2)]
					      + d3*(*field).u_comp[IND_2(k,j+3)]
					      + d4*(*field).u_comp[IND_2(k,j+4)]
					      )*h_2;
	  }
	  }


	  //-----------------------> AT AND CLOSE TO RIGHT BOUNDARY
	  
	  /* Value at n-1 */
	  (*dfield).u_comp[IND_2(k,n2f-1)] = -(d00*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d01*(*field).u_comp[IND_2(k,n2f-2)]
					      + d02*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d03*(*field).u_comp[IND_2(k,n2f-4)]
					      + d04*(*field).u_comp[IND_2(k,n2f-5)]
					      + d05*(*field).u_comp[IND_2(k,n2f-6)]
					      + d06*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d07*(*field).u_comp[IND_2(k,n2f-8)]
					      + d08*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d09*(*field).u_comp[IND_2(k,n2f-10)]
					      + d010*(*field).u_comp[IND_2(k,n2f-11)]
					      + d011*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;

	  /* Value at n-2 */
	  (*dfield).u_comp[IND_2(k,n2f-2)] =  -(d10*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d11*(*field).u_comp[IND_2(k,n2f-2)]
					      + d12*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d13*(*field).u_comp[IND_2(k,n2f-4)]
					      + d14*(*field).u_comp[IND_2(k,n2f-5)]
					      + d15*(*field).u_comp[IND_2(k,n2f-6)]
					      + d16*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d17*(*field).u_comp[IND_2(k,n2f-8)]
					      + d18*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d19*(*field).u_comp[IND_2(k,n2f-10)]
					      + d110*(*field).u_comp[IND_2(k,n2f-11)]
					      + d111*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;

	  /* Value at n-3 */
	  (*dfield).u_comp[IND_2(k,n2f-3)] =  -(d20*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d21*(*field).u_comp[IND_2(k,n2f-2)]
					      + d22*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d23*(*field).u_comp[IND_2(k,n2f-4)]
					      + d24*(*field).u_comp[IND_2(k,n2f-5)]
					      + d25*(*field).u_comp[IND_2(k,n2f-6)]
					      + d26*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d27*(*field).u_comp[IND_2(k,n2f-8)]
					      + d28*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d29*(*field).u_comp[IND_2(k,n2f-10)]
					      + d210*(*field).u_comp[IND_2(k,n2f-11)]
					      + d211*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;

	  /* Value at n-4 */
	  (*dfield).u_comp[IND_2(k,n2f-4)] = -(d30*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d31*(*field).u_comp[IND_2(k,n2f-2)]
					      + d32*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d33*(*field).u_comp[IND_2(k,n2f-4)]
					      + d34*(*field).u_comp[IND_2(k,n2f-5)]
					      + d35*(*field).u_comp[IND_2(k,n2f-6)]
					      + d36*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d37*(*field).u_comp[IND_2(k,n2f-8)]
					      + d38*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d39*(*field).u_comp[IND_2(k,n2f-10)]
					      + d310*(*field).u_comp[IND_2(k,n2f-11)]
					      + d311*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;

	  /* Value at n-5 */
	  (*dfield).u_comp[IND_2(k,n2f-5)] = -(d40*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d41*(*field).u_comp[IND_2(k,n2f-2)]
					      + d42*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d43*(*field).u_comp[IND_2(k,n2f-4)]
					      + d44*(*field).u_comp[IND_2(k,n2f-5)]
					      + d45*(*field).u_comp[IND_2(k,n2f-6)]
					      + d46*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d47*(*field).u_comp[IND_2(k,n2f-8)]
					      + d48*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d49*(*field).u_comp[IND_2(k,n2f-10)]
					      + d410*(*field).u_comp[IND_2(k,n2f-11)]
					      + d411*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;

	  /* Value at n-6 */
	  (*dfield).u_comp[IND_2(k,n2f-6)] = -(d50*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d51*(*field).u_comp[IND_2(k,n2f-2)]
					      + d52*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d53*(*field).u_comp[IND_2(k,n2f-4)]
					      + d54*(*field).u_comp[IND_2(k,n2f-5)]
					      + d55*(*field).u_comp[IND_2(k,n2f-6)]
					      + d56*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d57*(*field).u_comp[IND_2(k,n2f-8)]
					      + d58*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d59*(*field).u_comp[IND_2(k,n2f-10)]
					      + d510*(*field).u_comp[IND_2(k,n2f-11)]
					      + d511*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;


	  /* Value at n-7 */
	  (*dfield).u_comp[IND_2(k,n2f-7)] =  -(d60*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d61*(*field).u_comp[IND_2(k,n2f-2)]
					      + d62*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d63*(*field).u_comp[IND_2(k,n2f-4)]
					      + d64*(*field).u_comp[IND_2(k,n2f-5)]
					      + d65*(*field).u_comp[IND_2(k,n2f-6)]
					      + d66*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d67*(*field).u_comp[IND_2(k,n2f-8)]
					      + d68*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d69*(*field).u_comp[IND_2(k,n2f-10)]
					      + d610*(*field).u_comp[IND_2(k,n2f-11)]
					      + d611*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;


	  /* Value at n-8 */
	  (*dfield).u_comp[IND_2(k,n2f-8)] =  -(d70*(*field).u_comp[IND_2(k,n2f-1)] 
					      + d71*(*field).u_comp[IND_2(k,n2f-2)]
					      + d72*(*field).u_comp[IND_2(k,n2f-3)] 
					      + d73*(*field).u_comp[IND_2(k,n2f-4)]
					      + d74*(*field).u_comp[IND_2(k,n2f-5)]
					      + d75*(*field).u_comp[IND_2(k,n2f-6)]
					      + d76*(*field).u_comp[IND_2(k,n2f-7)] 
					      + d77*(*field).u_comp[IND_2(k,n2f-8)]
					      + d78*(*field).u_comp[IND_2(k,n2f-9)] 
					      + d79*(*field).u_comp[IND_2(k,n2f-10)]
					      + d710*(*field).u_comp[IND_2(k,n2f-11)]
					      + d711*(*field).u_comp[IND_2(k,n2f-12)]
					      )*h_2;

	  
	
      }
      }


#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
#ifdef DEBUG
      printf("at case 2");
#endif
    }

	
      break;

	
    default:{
      printf("wrong value given for coordinate\n");
      exit(0);
    }
      }
  //  printf("h_1= %f\n n1i= %d\n n1f= %d\n",h_1,n1i,n1f);
  //  printf("h_2= %f\n n2i= %d\n n2f= %d\n",h_2,n2i,n2f);
}



/* ------------------------------------------------------------------------------------------ */




/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */

#ifdef DISSIPATION

/* Second order dissipation operator with u_{-1} = u_{0}                                      */

void diss_S_2d(struct GRID_PAR *grid,
	     struct field *field,
	     struct field *dfield)
{/* diss_S_2d */

    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;
    int n_gridpts_1 = n1f-n1i;
    int n_gridpts_2 = n2f-n2i;

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;

    FLOAT d2 =  2.0;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);



	/* inner points */
	{register int ind1;
	register int ind2;
	  for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

		(*dfield).u_comp[IND_2(ind1,ind2)] = ((*field).u_comp[IND_2(ind1-1,ind2)]
					    	 - d2*(*field).u_comp[IND_2(ind1,ind2)]
					    	   +  (*field).u_comp[IND_2(ind1+1,ind2)]
						      )*h_1
					           + ((*field).u_comp[IND_2(ind1,ind2-1)]
					    	 - d2*(*field).u_comp[IND_2(ind1,ind2)]
					    	   +  (*field).u_comp[IND_2(ind1,ind2+1)]
						      )*h_2;
		}
		}
	}

 	/* x=0,1 faces */

	{register int ind2;
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){
		(*dfield).u_comp[IND_2(n1i,ind2)] = -d2*((*field).u_comp[IND_2(n1i,ind2)]
					    	       - (*field).u_comp[IND_2(n1i+1,ind2)]
							 )*h_1

						+ ((*field).u_comp[IND_2(n1i,ind2-1)]
					      - d2*(*field).u_comp[IND_2(n1i,ind2)]
					    	+  (*field).u_comp[IND_2(n1i,ind2+1)]
						   )*h_2;


		(*dfield).u_comp[IND_2(n1f-1,ind2)] = -d2*((*field).u_comp[IND_2(n1f-1,ind2)]
					    		 - (*field).u_comp[IND_2(n1f-2,ind2)]
							   )*h_1

						+ ((*field).u_comp[IND_2(n1f-1,ind2-1)]
					      - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
					       	+  (*field).u_comp[IND_2(n1f-1,ind2+1)]
						   )*h_2;
		}
	}

	 	/* y=0,1 faces */

	{register int ind1;
	    for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
		(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
					    	- d2*(*field).u_comp[IND_2(ind1,n2i)]
					          +  (*field).u_comp[IND_2(ind1+1,n2i)]
						     )*h_1

						- d2*((*field).u_comp[IND_2(ind1,n2i)]
					    	    - (*field).u_comp[IND_2(ind1,n2i+1)]
						      )*h_2;



		(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
					    	  - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
					    	    +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
						       )*h_1

						 - d2*((*field).u_comp[IND_2(ind1,n2f-1)]
					    	     - (*field).u_comp[IND_2(ind1,n2f-2)]
					               )*h_2;
		}
	}


 /* edges */
 
 /* runs on x */

	{register int ind1;
	    for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
			       - d2*(*field).u_comp[IND_2(ind1,n2i)]
				 +  (*field).u_comp[IND_2(ind1+1,n2i)]
				    )*h_1

			      + d2*((*field).u_comp[IND_2(ind1,n2i+1)]
				  - (*field).u_comp[IND_2(ind1,n2i)]
				    )*h_2;


(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
				 - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
				   +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
				      )*h_1

				- d2*((*field).u_comp[IND_2(ind1,n2f-1)]
				    - (*field).u_comp[IND_2(ind1,n2f-2)]
				      )*h_2;
	    }
	}



/* runs on y */

	{register int ind2;
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

(*dfield).u_comp[IND_2(n1i,ind2)] = d2*(
				       (*field).u_comp[IND_2(n1i+1,ind2)]
				     - (*field).u_comp[IND_2(n1i,ind2)]
				       )*h_1

				    + (
				       (*field).u_comp[IND_2(n1i,ind2-1)]
				  - d2*(*field).u_comp[IND_2(n1i,ind2)]
				     + (*field).u_comp[IND_2(n1i,ind2+1)]
				       )*h_2;


(*dfield).u_comp[IND_2(n1f-1,ind2)] = - d2*(
					   (*field).u_comp[IND_2(n1f-1,ind2)]
					 - (*field).u_comp[IND_2(n1f-2,ind2)]
					   )*h_1

				        + (
					   (*field).u_comp[IND_2(n1f-1,ind2-1)]
				      - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
					 + (*field).u_comp[IND_2(n1f-1,ind2+1)]
					   )*h_2;
	    }
	}




/* corners */

(*dfield).u_comp[IND_2(n1i,n2i)] = d2*((*field).u_comp[IND_2(n1i+1,n2i)]
				    - (*field).u_comp[IND_2(n1i,n2i)])*h_1

			        + d2*((*field).u_comp[IND_2(n1i,n2i+1)]
				    - (*field).u_comp[IND_2(n1i,n2i)])*h_2;

(*dfield).u_comp[IND_2(n1i,n2f-1)] = d2*((*field).u_comp[IND_2(n1i+1,n2f-1)]
				      - (*field).u_comp[IND_2(n1i,n2f-1)])*h_1

				  + d2*((*field).u_comp[IND_2(n1i,n2f-2)]
				      - (*field).u_comp[IND_2(n1i,n2f-1)])*h_2;


(*dfield).u_comp[IND_2(n1f-1,n2i)] = d2*((*field).u_comp[IND_2(n1f-2,n2i)]
				      - (*field).u_comp[IND_2(n1f-1,n2i)])*h_1

				  + d2*((*field).u_comp[IND_2(n1f-1,n2i+1)]
				      - (*field).u_comp[IND_2(n1f-1,n2i)])*h_2;

(*dfield).u_comp[IND_2(n1f-1,n2f-1)] = d2*((*field).u_comp[IND_2(n1f-2,n2f-1)]
				        - (*field).u_comp[IND_2(n1f-1,n2f-1)])*h_1

				    + d2*((*field).u_comp[IND_2(n1f-1,n2f-2)]
				        - (*field).u_comp[IND_2(n1f-1,n2f-1)])*h_2;


}


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */





/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */

/* Second order dissipation operator with u_{-1} = 2u_{0} - u_{1}                             */
/* This is just cero on the boundaries */

void diss_S0_2d(struct GRID_PAR *grid,
	     struct field *field,
	     struct field *dfield)
{

    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;

    int n_gridpts_1 = n1f-n1i;
    int n_gridpts_2 = n2f-n2i;

#ifdef EXCISION
    int exc_n1 = (*grid).exc_grid_1;
    int exc_n2 = (*grid).exc_grid_2;
    int exc_n1i = (n1f-n1i-1)/2 - exc_n1; 
    int exc_n2i = (n2f-n2i-1)/2 - exc_n2;
    int exc_n1f = (n1f-n1i-1)/2 + exc_n1; 
    int exc_n2f = (n2f-n2i-1)/2 + exc_n2;
#endif

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;


    FLOAT d2 =  2.0;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);



	/* inner points */
	{register int ind1;
	register int ind2;
	  for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

		(*dfield).u_comp[IND_2(ind1,ind2)] = ((*field).u_comp[IND_2(ind1-1,ind2)]
					    	 - d2*(*field).u_comp[IND_2(ind1,ind2)]
					    	   +  (*field).u_comp[IND_2(ind1+1,ind2)]
						      )*h_1
					           + ((*field).u_comp[IND_2(ind1,ind2-1)]
					    	 - d2*(*field).u_comp[IND_2(ind1,ind2)]
					    	   +  (*field).u_comp[IND_2(ind1,ind2+1)]
						      )*h_2
					          ;
		}
		}
	}


#ifdef EXCISION
	/* set to zero all second derivatives at interior points */
	  {register int n1;
	  register int n2;
	  for (n1=exc_n1i+1; n1 < exc_n1f; n1++){
	    for (n2=exc_n2i+1; n2 < exc_n2f; n2++){
		(*dfield).u_comp[IND_2(n1,n2)] = 0.0;
	    }
	  }
	  }
#endif	 


 	/* x=0,1 faces */

	{register int ind2;
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

		(*dfield).u_comp[IND_2(n1i,ind2)] =  ((*field).u_comp[IND_2(n1i,ind2-1)]
					         - d2*(*field).u_comp[IND_2(n1i,ind2)]
					    	   +  (*field).u_comp[IND_2(n1i,ind2+1)]
						      )*h_2;


		(*dfield).u_comp[IND_2(n1f-1,ind2)] =((*field).u_comp[IND_2(n1f-1,ind2-1)]
					         - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
					       	   +  (*field).u_comp[IND_2(n1f-1,ind2+1)]
						      )*h_2;


		}
	}


#ifdef EXCISION
	/* set to zero all second derivatives at interior points */
	  {register int ind2;
	    for (ind2 = exc_n2i+1; ind2 < exc_n2f; ind2++){

		(*dfield).u_comp[IND_2(exc_n1i,ind2)] =  ((*field).u_comp[IND_2(exc_n1i,ind2)]
					         - d2*(*field).u_comp[IND_2(exc_n1i,ind2)]
					    	   +  (*field).u_comp[IND_2(exc_n1i,ind2+1)]
						      )*h_2;


		(*dfield).u_comp[IND_2(exc_n1f,ind2)] =((*field).u_comp[IND_2(exc_n1f,ind2-1)]
					         - d2*(*field).u_comp[IND_2(exc_n1f,ind2)]
					       	   +  (*field).u_comp[IND_2(exc_n1f,ind2+1)]
						      )*h_2;
	    }
	  }
	  
#endif	 




	 	/* y=0,1 faces */

	{register int ind1;
	    for (ind1 = n1i+1; ind1 < n1f-1; ind1++){

		(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
					    	- d2*(*field).u_comp[IND_2(ind1,n2i)]
					          +  (*field).u_comp[IND_2(ind1+1,n2i)]
						     )*h_1;



		(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
					    	  - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
					    	    +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
						       )*h_1;

		}
	}

#ifdef EXCISION
	/* set to zero all second derivatives at interior points */
	  {register int ind1;
	    for (ind1 = exc_n1i+1; ind1 < exc_n1f; ind1++){

		(*dfield).u_comp[IND_2(ind1,exc_n2i)] = ((*field).u_comp[IND_2(ind1-1,exc_n2i)]
					    	- d2*(*field).u_comp[IND_2(ind1,exc_n2i)]
					          +  (*field).u_comp[IND_2(ind1+1,exc_n2i)]
						     )*h_1;



		(*dfield).u_comp[IND_2(ind1,exc_n2f)] = ((*field).u_comp[IND_2(ind1-1,exc_n2f)]
					    	  - d2*(*field).u_comp[IND_2(ind1,exc_n2f)]
					    	    +  (*field).u_comp[IND_2(ind1+1,exc_n2f)]
						       )*h_1;
	    }
	  }
	  
#endif	 






 /* edges */
 
 /* runs on x */

	{register int ind1;
	    for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
			       - d2*(*field).u_comp[IND_2(ind1,n2i)]
				 +  (*field).u_comp[IND_2(ind1+1,n2i)]
				    )*h_1;


(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
				 - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
				   +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
				      )*h_1;

				


(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
				 - d2*(*field).u_comp[IND_2(ind1,n2i)]
				   +  (*field).u_comp[IND_2(ind1+1,n2i)]
				      )*h_1;



(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
				   - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
				     +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
					)*h_1;


						}
	}



/* runs on y */

	{register int ind2;
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

(*dfield).u_comp[IND_2(n1i,ind2)] = (
				    (*field).u_comp[IND_2(n1i,ind2-1)]
			       - d2*(*field).u_comp[IND_2(n1i,ind2)]
			       +    (*field).u_comp[IND_2(n1i,ind2+1)]
				    )*h_2;


(*dfield).u_comp[IND_2(n1f-1,ind2)] = (
				    (*field).u_comp[IND_2(n1f-1,ind2-1)]
			       - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
			       +    (*field).u_comp[IND_2(n1f-1,ind2+1)]
				    )*h_2;



(*dfield).u_comp[IND_2(n1i,ind2)] = (
				    (*field).u_comp[IND_2(n1i,ind2-1)]
			       - d2*(*field).u_comp[IND_2(n1i,ind2)]
			       +    (*field).u_comp[IND_2(n1i,ind2+1)]
				    )*h_2;


(*dfield).u_comp[IND_2(n1f-1,ind2)] =(
				    (*field).u_comp[IND_2(n1f-1,ind2-1)]
			       - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
			       +    (*field).u_comp[IND_2(n1f-1,ind2+1)]
				    )*h_2;

						}
	}


	

/* corners */

(*dfield).u_comp[IND_2(n1i,n2i)] = 0.0;

(*dfield).u_comp[IND_2(n1i,n2i)] = 0.0;

(*dfield).u_comp[IND_2(n1i,n2f-1)] = 0.0;

(*dfield).u_comp[IND_2(n1i,n2f-1)] = 0.0;

(*dfield).u_comp[IND_2(n1f-1,n2i)] = 0.0;

(*dfield).u_comp[IND_2(n1f-1,n2i)] = 0.0;

(*dfield).u_comp[IND_2(n1f-1,n2f-1)] = 0.0;

(*dfield).u_comp[IND_2(n1f-1,n2f-1)] = 0.0;

}


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------------------------ */


/* ----------------------------------- DISSIPATIVE OPERATORS -------------------------------- */

/* ------------------------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------------------------ */




//-------------------------- KO-type dissipation, fourth derivative in the interior --------------------------//
// Written by Manuel.
void diss_KO4_2d(struct GRID_PAR *grid,
	      struct field *field,
		 struct field *dfield)
{/* diss_KO4 */

	 
 

  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;


  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  FLOAT d2 = 2.0;
  FLOAT d4 = 4.0;
  FLOAT d5 = 5.0;
  FLOAT d6 = 6.0;

  FLOAT p1;
  FLOAT p2;

  FLOAT h_1 = (xf-xi)/((FLOAT)(n1f-n1i-1)); // this is dx
  FLOAT h_2 = (yf-yi)/((FLOAT)(n2f-n2i-1)); // dy

  int n1r;
  int n1s;
  int n2r;
  int n2s;
  
  // on the sphere
  FLOAT sigma_sphere_rm2; // scalar product at gridpoint R Minus 2
  FLOAT sigma_sphere_rm1; // scalar product at gridpoint R Minus 1
  FLOAT sigma_sphere_sp2; // scalar product at gridpoint S Plus 2
  FLOAT sigma_sphere_sp1; // scalar product at gridpoint S Plus 1
  


  // define what scalar product, interior means (ie which are gridpoints r,s), depending on the derivatives that we are using

  char macro_value_strg[100];
  GET_MACRO_VALUE(DERIV);

  // --------- derivatives on the sphere ------------
  if (strcmp(macro_value_strg,"derivS_2d")==0) { // second order in the interior
    n1r = n1i+2;  n1s = n1f-2;  n2r = n2i+2;  n2s = n2f-2;
    p1 = 1.0/h_1;  p2 = 1.0/h_2; 
    sigma_sphere_rm2 = 0.5;  sigma_sphere_rm1 = 1.0;
  }
  else if (strcmp(macro_value_strg,"derivQ_2d")==0) {// fourth order in the interior
    n1r = n1i+4;  n1s = n1f-4;  n2r = n2i+4;  n2s = n2f-4;
    p1 = h_1;  p2 = h_2;
    sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;  
  }
   else if (strcmp(macro_value_strg,"derivQ_Per_x_2d")==0) {// fourth order in the interior
    n1r = n1i+4;  n1s = n1f-4;  n2r = n2i+4;  n2s = n2f-4;
    p1 = h_1;  p2 = h_2;
    sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;  
  } else if (strcmp(macro_value_strg,"derivQ_Per_y_2d")==0) {// fourth order in the interior
    n1r = n1i+4;  n1s = n1f-4;  n2r = n2i+4;  n2s = n2f-4;
    p1 = h_1;  p2 = h_2;
    sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;  
  }
  else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_2d")==0) { // sixth order in the interior
    n1r = n1i+6;  n1s = n1f-6;  n2r = n2i+6;  n2s = n2f-6;
    p1 = pow(h_1,3); p2 = pow(h_2,3);
    sigma_sphere_rm2 = 7877.0/8640.0;  sigma_sphere_rm1 = 43801.0/43200.0;  
  }
  else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_2d")==0) {// eighth order in the interior
  n1r = n1i+8;  n1s = n1f-8;  n2r = n2i+8;  n2s = n2f-8;
  p1 = pow(h_1,5);  p2 = pow(h_2,5);
  sigma_sphere_rm2 = 670091.0/725760.0;  sigma_sphere_rm1 = 5127739.0/5080320.0;  
  }
  else {
    printf("I don't know what derivative you are using, will stop here.\n");
    exit(1);
  }
  
  sigma_sphere_sp2 = sigma_sphere_rm2;  sigma_sphere_sp1 = sigma_sphere_rm1;


#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif


  // initialize to zero
  {register int ind1;
  register int ind2;
  register int ind3;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	(*dfield).u_comp[IND_2(ind1,ind2)] =  0.0;
      }
    }
  }
  
  
 
  //------------------------- dissipation in the x-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1r; ind1 < n1s; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  -((*field).u_comp[IND_2(ind1-2,ind2)]  - 
						     d4*(*field).u_comp[IND_2(ind1-1,ind2)]  + 
						     d6*(*field).u_comp[IND_2(ind1+0,ind2)] - 
						     d4*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
						     (*field).u_comp[IND_2(ind1+2,ind2)])*p1;
      }
    }
  }
  
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind2 = n2i; ind2 < n2f; ind2++){
      
      // left one
      ind1 = n1r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  -((*field).u_comp[IND_2(ind1+0,ind2)]  - 
						      d2*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
						      (*field).u_comp[IND_2(ind1+2,ind2)])*p1/sigma_sphere_rm2;
      
      ind1 = n1r-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  -( - d2*(*field).u_comp[IND_2(ind1-1,ind2)]  + 
						    d5*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
						    d4*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
						    (*field).u_comp[IND_2(ind1+2,ind2)])*p1/sigma_sphere_rm1;
      
      
      // right one
      ind1 = n1s;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  -( - d2*(*field).u_comp[IND_2(ind1+1,ind2)]  + 
						    d5*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
						    d4*(*field).u_comp[IND_2(ind1-1,ind2)]  +    
						    (*field).u_comp[IND_2(ind1-2,ind2)])*p1/sigma_sphere_sp1;
      
      ind1 = n1s+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] = -((*field).u_comp[IND_2(ind1+0,ind2)]  - 
						     d2*(*field).u_comp[IND_2(ind1-1,ind2)]  +  
						     (*field).u_comp[IND_2(ind1-2,ind2)])*p1/sigma_sphere_sp2;
      
      
    
  }
  }
  
  //------------------------- add dissipation in the y-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2r; ind2 < n2s; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)] 
	  -((*field).u_comp[IND_2(ind1,ind2-2)]  - 
	    d4*(*field).u_comp[IND_2(ind1,ind2-1)]  + 
	    d6*(*field).u_comp[IND_2(ind1,ind2+0)] - 
	    d4*(*field).u_comp[IND_2(ind1,ind2+1)]  +  
	    (*field).u_comp[IND_2(ind1,ind2+2)])*p2;
      
    }
  }
  }
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
      
      // left one
      ind2 = n2r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	-((*field).u_comp[IND_2(ind1,ind2+0)]  - 
	     d2*(*field).u_comp[IND_2(ind1,ind2+1)]  +  
	     (*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_rm2 ;
      
      ind2 = n2r-1; 
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	-( - d2*(*field).u_comp[IND_2(ind1,ind2-1)]  + 
	   d5*(*field).u_comp[IND_2(ind1,ind2+0)]  - 
	   d4*(*field).u_comp[IND_2(ind1,ind2+1)]  +  
	   (*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_rm1 ;
      
      
      // right one
      ind2 = n2s; 
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	-( - d2*(*field).u_comp[IND_2(ind1,ind2+1)]  + 
	   d5*(*field).u_comp[IND_2(ind1,ind2+0)]  - 
	   d4*(*field).u_comp[IND_2(ind1,ind2-1)]  +    
	   (*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_sp1 ;
      
      ind2 = n2s+1; 
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	-((*field).u_comp[IND_2(ind1,ind2+0)]  - 
	     d2*(*field).u_comp[IND_2(ind1,ind2-1)]  +  
	     (*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_sp2 ;
      
      
    }
  }
  
  
  

  
}



//-------------------------- KO-type dissipation, fourth derivative in the interior --------------------------//
// Written by Manuel.
void diss_KO4_Per_x_2d(struct GRID_PAR *grid,
	      struct field *field,
		 struct field *dfield)
{/* diss_KO4 */

	 
 

  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;


  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  FLOAT d2 = 2.0;
  FLOAT d4 = 4.0;
  FLOAT d5 = 5.0;
  FLOAT d6 = 6.0;

  FLOAT p1;
  FLOAT p2;

  FLOAT h_1 = (xf-xi)/((FLOAT)(n1f-n1i)); // this is dx
  FLOAT h_2 = (yf-yi)/((FLOAT)(n2f-n2i-1)); // dy

  int n1r;
  int n1s;
  int n2r;
  int n2s;
  
  // on the sphere
  FLOAT sigma_sphere_rm2; // scalar product at gridpoint R Minus 2
  FLOAT sigma_sphere_rm1; // scalar product at gridpoint R Minus 1
  FLOAT sigma_sphere_sp2; // scalar product at gridpoint S Plus 2
  FLOAT sigma_sphere_sp1; // scalar product at gridpoint S Plus 1
  


  // define what scalar product, interior means (ie which are gridpoints r,s), depending on the derivatives that we are using

  char macro_value_strg[100];
  GET_MACRO_VALUE(DERIV);

  // --------- derivatives on the sphere ------------
  if (strcmp(macro_value_strg,"derivS_Per_x_2d")==0) { // second order in the interior
    n1r = n1i+2;  n1s = n1f-2;  n2r = n2i+2;  n2s = n2f-2;
    p1 = 1.0/h_1;  p2 = 1.0/h_2; 
    sigma_sphere_rm2 = 0.5;  sigma_sphere_rm1 = 1.0;

  }
  else if (strcmp(macro_value_strg,"derivS_Per_x_2d")==0) { // second order in the interior
    n1r = n1i+2;  n1s = n1f-2;  n2r = n2i+2;  n2s = n2f-2;
    p1 = 1.0/h_1;  p2 = 1.0/h_2; 
    sigma_sphere_rm2 = 0.5;  sigma_sphere_rm1 = 1.0;
  }
   else if (strcmp(macro_value_strg,"derivQ_Per_x_2d")==0) {// fourth order in the interior
    n1r = n1i+4;  n1s = n1f-4;  n2r = n2i+4;  n2s = n2f-4;
    p1 = h_1;  p2 = h_2;
    sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;  
  }

  else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_Per_x_2d")==0) { // sixth order in the interior
    n1r = n1i+6;  n1s = n1f-6;  n2r = n2i+6;  n2s = n2f-6;
    p1 = pow(h_1,3); p2 = pow(h_2,3);
    sigma_sphere_rm2 = 7877.0/8640.0;  sigma_sphere_rm1 = 43801.0/43200.0;  
  }
  else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_Per_x_2d")==0) {// eighth order in the interior
  n1r = n1i+8;  n1s = n1f-8;  n2r = n2i+8;  n2s = n2f-8;
  p1 = pow(h_1,5);  p2 = pow(h_2,5);
  sigma_sphere_rm2 = 670091.0/725760.0;  sigma_sphere_rm1 = 5127739.0/5080320.0;  
  }
  else {
    printf("I don't know what derivative you are using, will stop here. Look at derivs_2d.c \n");
    exit(1);
  }
  
  sigma_sphere_sp2 = sigma_sphere_rm2;  sigma_sphere_sp1 = sigma_sphere_rm1;


#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif


  // initialize to zero
  {register int ind1;
  register int ind2;
  register int ind3;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	(*dfield).u_comp[IND_2(ind1,ind2)] =  0.0;
      }
    }
  }
  
  
 
  //------------------------- dissipation in the x-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1i+2; ind1 < n1f-2; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  -((*field).u_comp[IND_2(ind1-2,ind2)]  - 
				             d4*(*field).u_comp[IND_2(ind1-1,ind2)]  + 
				             d6*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
				             d4*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
					        (*field).u_comp[IND_2(ind1+2,ind2)])*p1;
      }
    }
  }
  
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind2 = n2i; ind2 < n2f; ind2++){
      
      // left one
      ind1 = n1i;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  -((*field).u_comp[IND_2(n1f -2,ind2)]  - 
				           d4*(*field).u_comp[IND_2(n1f -1,ind2)]  + 
				           d6*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
				           d4*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
					      (*field).u_comp[IND_2(ind1+2,ind2)])*p1;
      
      ind1 = n1i+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  -((*field).u_comp[IND_2(n1f -1,ind2)]  - 
					   d4*(*field).u_comp[IND_2(ind1-1,ind2)]  + 
				           d6*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
					   d4*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
					      (*field).u_comp[IND_2(ind1+2,ind2)])*p1;
      
      
      // right one
      ind1 = n1f-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  -((*field).u_comp[IND_2(ind1-2,ind2)]  - 
					   d4*(*field).u_comp[IND_2(ind1-1,ind2)]  + 
				           d6*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
				           d4*(*field).u_comp[IND_2(n1i   ,ind2)]  +  
				       	      (*field).u_comp[IND_2(n1i +1,ind2)])*p1;
      
      
      ind1 = n1f-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = -((*field).u_comp[IND_2(ind1-2,ind2)]  - 
				          d4*(*field).u_comp[IND_2(ind1-1,ind2)]  + 
			                  d6*(*field).u_comp[IND_2(ind1+0,ind2)]  - 
				          d4*(*field).u_comp[IND_2(ind1+1,ind2)]  +  
				             (*field).u_comp[IND_2(n1i   ,ind2)])*p1;
      
      
    
  }
  }
  
  //------------------------- add dissipation in the y-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2r; ind2 < n2s; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)] 
	  -((*field).u_comp[IND_2(ind1,ind2-2)]  - 
	    d4*(*field).u_comp[IND_2(ind1,ind2-1)]  + 
	    d6*(*field).u_comp[IND_2(ind1,ind2+0)] - 
	    d4*(*field).u_comp[IND_2(ind1,ind2+1)]  +  
	    (*field).u_comp[IND_2(ind1,ind2+2)])*p2;
      
    }
  }
  }
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
      
      // left one
      ind2 = n2r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	-((*field).u_comp[IND_2(ind1,ind2+0)]  - 
	     d2*(*field).u_comp[IND_2(ind1,ind2+1)]  +  
	     (*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_rm2 ;
      
      ind2 = n2r-1; 
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	-( - d2*(*field).u_comp[IND_2(ind1,ind2-1)]  + 
	   d5*(*field).u_comp[IND_2(ind1,ind2+0)]  - 
	   d4*(*field).u_comp[IND_2(ind1,ind2+1)]  +  
	   (*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_rm1 ;
      
      
      // right one
      ind2 = n2s; 
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	-( - d2*(*field).u_comp[IND_2(ind1,ind2+1)]  + 
	   d5*(*field).u_comp[IND_2(ind1,ind2+0)]  - 
	   d4*(*field).u_comp[IND_2(ind1,ind2-1)]  +    
	   (*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_sp1 ;
      
      ind2 = n2s+1; 
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	-((*field).u_comp[IND_2(ind1,ind2+0)]  - 
	     d2*(*field).u_comp[IND_2(ind1,ind2-1)]  +  
	     (*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_sp2 ;
      
      
    }
  }
  
  
  

  
}


//-------------------------- KO-type dissipation, sixth derivative in the interior --------------------------//
// Written by Manuel.

void diss_KO6_2d(struct GRID_PAR *grid,
	      struct field *field,
		 struct field *dfield)
{/* diss_KO6 */



  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;


  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  FLOAT d2 = 2.0;
  FLOAT d4 = 4.0;
  FLOAT d5 = 5.0;
  FLOAT d6 = 6.0;

  FLOAT p1;
  FLOAT p2;

  FLOAT h_1 = (xf-xi)/((FLOAT)(n1f-n1i-1)); // this is dx
  FLOAT h_2 = (yf-yi)/((FLOAT)(n2f-n2i-1)); // dy

  int n1r;
  int n1s;
  int n2r;
  int n2s;

  // scalar product on the sphere
  FLOAT sigma_sphere_rm3; // scalar product at gridpoint R Minus 3  
  FLOAT sigma_sphere_rm2; // scalar product at gridpoint R Minus 2
  FLOAT sigma_sphere_rm1; // scalar product at gridpoint R Minus 1
  FLOAT sigma_sphere_sp3; // scalar product at gridpoint S Plus 3
  FLOAT sigma_sphere_sp2; // scalar product at gridpoint S Plus 2
  FLOAT sigma_sphere_sp1; // scalar product at gridpoint S Plus 1
 


  // define what scalar product, interior means (ie which are gridpoints r,s), depending on the derivatives that we are using

  char macro_value_strg[100];
  GET_MACRO_VALUE(DERIV);

  // --------- derivatives on the sphere ------------
  if (strcmp(macro_value_strg,"derivS_2d")==0) { // second order in the interior
    n1r = n1i+3;  n1s = n1f-4;  n2r = n2i+3;  n2s = n2f-4;
    p1 = 1.0/pow(h_1,3);  p2 = 1.0/pow(h_2,3); 
    sigma_sphere_rm3 = 0.5; sigma_sphere_rm2 = 1.0;  sigma_sphere_rm1 = 1.0;
  }

  else if (strcmp(macro_value_strg,"derivQ_2d")==0) {// fourth order in the interior
    n1r = n1i+4;  n1s = n1f-5;  n2r = n2i+4;  n2s = n2f-5;
    p1 = 1.0/h_1;  p2 = 1.0/h_2;
    sigma_sphere_rm3 = 59.0/48.0; sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;  
  }



  else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_2")==0) { // sixth order in the interior
    n1r = n1i+6;  n1s = n1f-7;  n2r = n2i+6;  n2s = n2f-7;
    p1 = h_1; p2 = h_2;
    sigma_sphere_rm3= 5359.0/4320.0; sigma_sphere_rm2 = 7877.0/8640.0;  sigma_sphere_rm1 = 43801.0/43200.0;  
  }

  else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_2")==0) {// eighth order in the interior
  n1r = n1i+8;  n1s = n1f-9;  n2r = n2i+8;  n2s = n2f-9;
  p1 = pow(h_1,3);  p2 = pow(h_2,3);
  sigma_sphere_rm3 = 103097.0/80640.0; sigma_sphere_rm2 = 670091.0/725760.0;  sigma_sphere_rm1 = 5127739.0/5080320.0;  
  }

  else {
    printf("I don't know what derivative you are using, will stop here.\n");
    exit(1);
  }
 
  sigma_sphere_sp3 = sigma_sphere_rm3; sigma_sphere_sp2 = sigma_sphere_rm2;  sigma_sphere_sp1 = sigma_sphere_rm1;



#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
  
  // initialize to zero
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  0.0;
	
      } } }
  
 
  //------------------------- dissipation in the x-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1r; ind1 < n1s+1; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(ind1+3,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     + (*field).u_comp[IND_2(ind1-3,ind2)])*p1;
      } } } 
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind2 = n2i; ind2 < n2f; ind2++){
      
      // left one
      ind1 = n1r-3;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ((*field).u_comp[IND_2(ind1+3,ind2)] 
						  - 3.0*(*field).u_comp[IND_2(ind1+2,ind2)]  
						  + 3.0*(*field).u_comp[IND_2(ind1+1,ind2)] 
						  - (*field).u_comp[IND_2(ind1+0,ind2)])*p1/sigma_sphere_rm3;
      
      ind1 = n1r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ((*field).u_comp[IND_2(ind1+3,ind2)] 
						  - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]
						  + 12.0*(*field).u_comp[IND_2(ind1+1,ind2)]
						  - 10.0*(*field).u_comp[IND_2(ind1+0,ind2)]
						  + 3.0*(*field).u_comp[IND_2(ind1-1,ind2)])*p1/sigma_sphere_rm2;
      
      ind1 = n1r-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ((*field).u_comp[IND_2(ind1+3,ind2)] 
						  - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]
						  + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]
						  - 19.0*(*field).u_comp[IND_2(ind1+0,ind2)]
						  + 12.0*(*field).u_comp[IND_2(ind1-1,ind2)]
						  - 3.0*(*field).u_comp[IND_2(ind1-2,ind2)])*p1/sigma_sphere_rm1;
      
      // right one
      ind1 = n1s+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] = ((*field).u_comp[IND_2(ind1-3,ind2)] 
						  - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						  + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]
						  - 19.0*(*field).u_comp[IND_2(ind1+0,ind2)]
						  + 12.0*(*field).u_comp[IND_2(ind1+1,ind2)]
						  - 3.0*(*field).u_comp[IND_2(ind1+2,ind2)])*p1/sigma_sphere_sp1;
      
      ind1 = n1s+2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = ((*field).u_comp[IND_2(ind1-3,ind2)] 
						  - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						  + 12.0*(*field).u_comp[IND_2(ind1-1,ind2)]
						  - 10.0*(*field).u_comp[IND_2(ind1+0,ind2)]
						  + 3.0*(*field).u_comp[IND_2(ind1+1,ind2)])*p1/sigma_sphere_sp2;


      ind1 = n1s+3;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ((*field).u_comp[IND_2(ind1-3,ind2)] 
						  - 3.0*(*field).u_comp[IND_2(ind1-2,ind2)]  
						  + 3.0*(*field).u_comp[IND_2(ind1-1,ind2)] 
						  - (*field).u_comp[IND_2(ind1+0,ind2)])*p1/sigma_sphere_sp3;
      
    } } 

  //------------------------- add dissipation in the y-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2r; ind2 < n2s+1; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] = 	(*dfield).u_comp[IND_2(ind1,ind2)] 
	  + ((*field).u_comp[IND_2(ind1,ind2+3)] 
	     - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]   
	     + 15.0*(*field).u_comp[IND_2(ind1,ind2+1)]   
	     - 20.0*(*field).u_comp[IND_2(ind1,ind2+0)]  
	     + 15.0*(*field).u_comp[IND_2(ind1,ind2-1)]   
	     - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	     + (*field).u_comp[IND_2(ind1,ind2-3)])*p2;
      
      } } } 
  
  // boundaries

  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
      
      // left one
      ind2 = n2r-3;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)] 
	+ ((*field).u_comp[IND_2(ind1,ind2+3)] 
	   - 3.0*(*field).u_comp[IND_2(ind1,ind2+2)]  
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	   - (*field).u_comp[IND_2(ind1,ind2+0)])*p2/sigma_sphere_rm3;
      
      ind2 = n2r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)] 
	+ ((*field).u_comp[IND_2(ind1,ind2+3)] 
	   - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]
	   + 12.0*(*field).u_comp[IND_2(ind1,ind2+1)]
	   - 10.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2-1)])*p2/sigma_sphere_rm2;
      
      ind2 = n2r-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	+((*field).u_comp[IND_2(ind1,ind2+3)] 
	  - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]
	  + 15.0*(*field).u_comp[IND_2(ind1,ind2+1)]
	  - 19.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	  + 12.0*(*field).u_comp[IND_2(ind1,ind2-1)]
	  - 3.0*(*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_rm1;
      
      // right one
      ind2 = n2s+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ((*field).u_comp[IND_2(ind1,ind2-3)] 
	   - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	   + 15.0*(*field).u_comp[IND_2(ind1,ind2-1)]
	   - 19.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	   + 12.0*(*field).u_comp[IND_2(ind1,ind2+1)]
	   - 3.0*(*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_sp1;
      
      ind2 = n2s+2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ((*field).u_comp[IND_2(ind1,ind2-3)] 
	   - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	   + 12.0*(*field).u_comp[IND_2(ind1,ind2-1)]
	   - 10.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2+1)])*p2/sigma_sphere_sp2;
      
      ind2 = n2s+3;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ((*field).u_comp[IND_2(ind1,ind2-3)] 
	   - 3.0*(*field).u_comp[IND_2(ind1,ind2-2)]  
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	   - (*field).u_comp[IND_2(ind1,ind2+0)])*p2/sigma_sphere_sp3;

    } } 
      
   
}




//-------------------------- KO-type dissipation, sixth derivative in the interior --------------------------//
// Written by Manuel.

void diss_KO6_Per_x_2d(struct GRID_PAR *grid,
	      struct field *field,
		 struct field *dfield)
{/* diss_KO6 */



  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;


  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  FLOAT d2 = 2.0;
  FLOAT d4 = 4.0;
  FLOAT d5 = 5.0;
  FLOAT d6 = 6.0;

  FLOAT p1;
  FLOAT p2;

  FLOAT h_1 = (xf-xi)/((FLOAT)(n1f-n1i)); // this is dx, periodic
  FLOAT h_2 = (yf-yi)/((FLOAT)(n2f-n2i-1)); // dy

  int n1r;
  int n1s;
  int n2r;
  int n2s;

  // scalar product on the sphere
  FLOAT sigma_sphere_rm3; // scalar product at gridpoint R Minus 3  
  FLOAT sigma_sphere_rm2; // scalar product at gridpoint R Minus 2
  FLOAT sigma_sphere_rm1; // scalar product at gridpoint R Minus 1
  FLOAT sigma_sphere_sp3; // scalar product at gridpoint S Plus 3
  FLOAT sigma_sphere_sp2; // scalar product at gridpoint S Plus 2
  FLOAT sigma_sphere_sp1; // scalar product at gridpoint S Plus 1
 


  // define what scalar product, interior means (ie which are gridpoints r,s), depending on the derivatives that we are using

  char macro_value_strg[100];
  GET_MACRO_VALUE(DERIV);

  // --------- derivatives on the sphere ------------
  if (strcmp(macro_value_strg,"derivS_Per_x_2d")==0) { // second order in the interior
    n1r = n1i+3;  n1s = n1f-4;  n2r = n2i+3;  n2s = n2f-4;
    p1 = 1.0/pow(h_1,3);  p2 = 1.0/pow(h_2,3); 
    sigma_sphere_rm3 = 0.5; sigma_sphere_rm2 = 1.0;  sigma_sphere_rm1 = 1.0;
  }

   else if (strcmp(macro_value_strg,"derivQ_Per_x_2d")==0) {// fourth order in the interior

    n1r = n1i+4;  n1s = n1f-5;  n2r = n2i+4;  n2s = n2f-5;
    p1 = 1.0/h_1;  p2 = 1.0/h_2;
    sigma_sphere_rm3 = 59.0/48.0; sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;  
  }

  else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_Per_x_2d")==0) { // sixth order in the interior
    n1r = n1i+6;  n1s = n1f-7;  n2r = n2i+6;  n2s = n2f-7;
    p1 = h_1; p2 = h_2;
    sigma_sphere_rm3= 5359.0/4320.0; sigma_sphere_rm2 = 7877.0/8640.0;  sigma_sphere_rm1 = 43801.0/43200.0;  
  }

  else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_Per_x_2d")==0) {// eighth order in the interior
  n1r = n1i+8;  n1s = n1f-9;  n2r = n2i+8;  n2s = n2f-9;
  p1 = pow(h_1,3);  p2 = pow(h_2,3);
  sigma_sphere_rm3 = 103097.0/80640.0; sigma_sphere_rm2 = 670091.0/725760.0;  sigma_sphere_rm1 = 5127739.0/5080320.0;  
  }

  else {
    printf("I don't know what derivative you are using, will stop here. Look at derivs_2d.c\n");
    exit(1);
  }
 
  sigma_sphere_sp3 = sigma_sphere_rm3; sigma_sphere_sp2 = sigma_sphere_rm2;  sigma_sphere_sp1 = sigma_sphere_rm1;



#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
  
  // initialize to zero
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  0.0;
	
      } } }
  
 
  //------------------------- dissipation in the x-direction ----------------------------//
  // interior

//printf("p1=%f, p2=%f \n",p1,p2);

  {register int ind1;
  register int ind2;
  for (ind1 = n1i+3; ind1 < n1f-4; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(ind1+3,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     + (*field).u_comp[IND_2(ind1-3,ind2)])*p1;
      } } } 
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind2 = n2i; ind2 < n2f; ind2++){
      
      // left one

	ind1 = n1i;
(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(ind1+3,ind2)] 
					- 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
			         	+ 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
					- 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
					+ 15.0*(*field).u_comp[IND_2(n1f-1,ind2)]   
					- 6.0*(*field).u_comp[IND_2(n1f-2,ind2)]
					+ (*field).u_comp[IND_2(n1f-3,ind2)])*p1;
       
	ind1 = n1i+1;
(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(ind1+3,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(n1f-1,ind2)]
						     + (*field).u_comp[IND_2(n1f-2,ind2)])*p1;
	ind1 = n1i+2;
(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(ind1+3,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     + (*field).u_comp[IND_2(n1f-1,ind2)])*p1;

	// right boundary


	ind1 = n1f-1;
(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(n1i+2,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(n1i+1,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(n1i,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     + (*field).u_comp[IND_2(ind1-3,ind2)])*p1;

	ind1 = n1f-2;
(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(n1i+1,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(n1i,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     + (*field).u_comp[IND_2(ind1-3,ind2)])*p1;

	ind1 = n1f-3;
(*dfield).u_comp[IND_2(ind1,ind2)] =  ( (*field).u_comp[IND_2(n1i,ind2)] 
						     - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
						     + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     - 20.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     + (*field).u_comp[IND_2(ind1-3,ind2)])*p1;


    } } 

  //------------------------- add dissipation in the y-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2r; ind2 < n2s+1; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] = 	(*dfield).u_comp[IND_2(ind1,ind2)] 
	  + ((*field).u_comp[IND_2(ind1,ind2+3)] 
	     - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]   
	     + 15.0*(*field).u_comp[IND_2(ind1,ind2+1)]   
	     - 20.0*(*field).u_comp[IND_2(ind1,ind2+0)]  
	     + 15.0*(*field).u_comp[IND_2(ind1,ind2-1)]   
	     - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	     + (*field).u_comp[IND_2(ind1,ind2-3)])*p2;
      
      } } } 
  
  // boundaries

  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
      
      // left one
      ind2 = n2r-3;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)] 
	+ ((*field).u_comp[IND_2(ind1,ind2+3)] 
	   - 3.0*(*field).u_comp[IND_2(ind1,ind2+2)]  
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	   - (*field).u_comp[IND_2(ind1,ind2+0)])*p2/sigma_sphere_rm3;
      
      ind2 = n2r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)] 
	+ ((*field).u_comp[IND_2(ind1,ind2+3)] 
	   - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]
	   + 12.0*(*field).u_comp[IND_2(ind1,ind2+1)]
	   - 10.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2-1)])*p2/sigma_sphere_rm2;
      
      ind2 = n2r-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	+((*field).u_comp[IND_2(ind1,ind2+3)] 
	  - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]
	  + 15.0*(*field).u_comp[IND_2(ind1,ind2+1)]
	  - 19.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	  + 12.0*(*field).u_comp[IND_2(ind1,ind2-1)]
	  - 3.0*(*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_rm1;
      
      // right one
      ind2 = n2s+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ((*field).u_comp[IND_2(ind1,ind2-3)] 
	   - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	   + 15.0*(*field).u_comp[IND_2(ind1,ind2-1)]
	   - 19.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	   + 12.0*(*field).u_comp[IND_2(ind1,ind2+1)]
	   - 3.0*(*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_sp1;
      
      ind2 = n2s+2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ((*field).u_comp[IND_2(ind1,ind2-3)] 
	   - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	   + 12.0*(*field).u_comp[IND_2(ind1,ind2-1)]
	   - 10.0*(*field).u_comp[IND_2(ind1,ind2+0)]
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2+1)])*p2/sigma_sphere_sp2;
      
      ind2 = n2s+3;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ((*field).u_comp[IND_2(ind1,ind2-3)] 
	   - 3.0*(*field).u_comp[IND_2(ind1,ind2-2)]  
	   + 3.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	   - (*field).u_comp[IND_2(ind1,ind2+0)])*p2/sigma_sphere_sp3;

    } } 
      
   
}


//-------------------------- KO-type dissipation, eight derivative in the interior --------------------------//
// Written by Manuel.

void diss_KO8_2d(struct GRID_PAR *grid,
	      struct field *field,
		 struct field *dfield)
{ /* diss_KO8 */

  int n_gridpts_1 = (*grid).n_gridpts_1;
  int n_gridpts_2 = (*grid).n_gridpts_2;


  int n1i = (*grid).start_grid_1;
  int n1f = (*grid).final_grid_1;
  int n2i = (*grid).start_grid_2;
  int n2f = (*grid).final_grid_2;

  FLOAT xi = (*grid).initial_x;
  FLOAT xf = (*grid).final_x;
  FLOAT yi = (*grid).initial_y;
  FLOAT yf = (*grid).final_y;

  FLOAT p1;
  FLOAT p2;

  FLOAT h_1 = (xf-xi)/((FLOAT)(n1f-n1i-1)); // this is dx
  FLOAT h_2 = (yf-yi)/((FLOAT)(n2f-n2i-1)); // dy

  int n1r;
  int n1s;
  int n2r;
  int n2s;

  // scalar product on the sphere
  FLOAT sigma_sphere_rm4; // scalar product at gridpoint R Minus 3  
  FLOAT sigma_sphere_rm3; // scalar product at gridpoint R Minus 3  
  FLOAT sigma_sphere_rm2; // scalar product at gridpoint R Minus 2
  FLOAT sigma_sphere_rm1; // scalar product at gridpoint R Minus 1

  FLOAT sigma_sphere_sp4; // scalar product at gridpoint S Plus 3 4
  FLOAT sigma_sphere_sp3; // scalar product at gridpoint S Plus 3
  FLOAT sigma_sphere_sp2; // scalar product at gridpoint S Plus 2
  FLOAT sigma_sphere_sp1; // scalar product at gridpoint S Plus 1
  
 


  // define what scalar product, interior means (ie which are gridpoints r,s), depending on the derivatives that we are using

  char macro_value_strg[100];
  GET_MACRO_VALUE(DERIV);

  // --------- derivatives on the sphere ------------
  if (strcmp(macro_value_strg,"derivS_2d")==0) { // second order in the interior

    n1r = n1i+4;  n1s = n1f-5;  n2r = n2i+4;  n2s = n2f-5;
    p1 = 1.0/pow(h_1,5);  p2 = 1.0/pow(h_2,5); 
    
    sigma_sphere_rm4 = 0.5; sigma_sphere_rm3 = 1.0; sigma_sphere_rm2 = 1.0;  sigma_sphere_rm1 = 1.0;
  }
  
  else if (strcmp(macro_value_strg,"derivQ_2d")==0) {// fourth order in the interior

    n1r = n1i+4;  n1s = n1f-5;  n2r = n2i+4;  n2s = n2f-5;
    p1 = 1.0/pow(h_1,3);  p2 = 1.0/pow(h_2,3);
    
    sigma_sphere_rm4 = 17.0/48.0; sigma_sphere_rm3 = 59.0/48.0; sigma_sphere_rm2 = 43.0/48.0;  sigma_sphere_rm1 = 49.0/48.0;
  }

  else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_2d")==0) { // sixth order in the interior

    n1r = n1i+6;  n1s = n1f-7;  n2r = n2i+6;  n2s = n2f-7;
    p1 = 1.0/h_1; p2 = 1.0/h_2;

    sigma_sphere_rm4 = 2711.0/4320.0; sigma_sphere_rm3 = 5359.0/4320.0; 
    sigma_sphere_rm2 = 7877.0/8640.0;  sigma_sphere_rm1 = 43801.0/43200.0;
  }

  else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_2d")==0) {// eighth order in the interior

  n1r = n1i+8;  n1s = n1f-9;  n2r = n2i+8;  n2s = n2f-9;
  p1 = h_1;  p2 = h_2;
  
  sigma_sphere_rm4 = 299527.0/725760.0; sigma_sphere_rm3 = 103097.0/80640.0; 
  sigma_sphere_rm2 = 670091.0/725760.0;  sigma_sphere_rm1 = 5127739.0/5080320.0;
  }

  else {
    printf("I don't know what derivative you are using, will stop here.\n");
    exit(1);
  }
 
  sigma_sphere_sp4 = sigma_sphere_rm4;  sigma_sphere_sp3 = sigma_sphere_rm3;
  sigma_sphere_sp2 = sigma_sphere_rm2;  sigma_sphere_sp1 = sigma_sphere_rm1;
 



#ifdef EXCISION       
  printf("************ THIS DIFFERENCE OPERATOR DOES NOT WORK WITH EXCISION !!! **************");
#endif	 
  
  // initialize to zero
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	
	(*dfield).u_comp[IND_2(ind1,ind2)] =  0.0;
	
      } } } 
  
 
  //------------------------- dissipation in the x-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1r; ind1 < n1s+1; ind1++){
    for (ind2 = n2i; ind2 < n2f; ind2++){
	(*dfield).u_comp[IND_2(ind1,ind2)] =  -( (*field).u_comp[IND_2(ind1+4,ind2)] 
						     - 8.0*(*field).u_comp[IND_2(ind1+3,ind2)] 
						     + 28.0*(*field).u_comp[IND_2(ind1+2,ind2)]   
						     - 56.0*(*field).u_comp[IND_2(ind1+1,ind2)]   
						     + 70.0*(*field).u_comp[IND_2(ind1+0,ind2)]  
						     - 56.0*(*field).u_comp[IND_2(ind1-1,ind2)]   
						     + 28.0*(*field).u_comp[IND_2(ind1-2,ind2)]
						     - 8.0*(*field).u_comp[IND_2(ind1-3,ind2)]
						     + (*field).u_comp[IND_2(ind1-4,ind2)]
						     )*p1;
      } } } 
  
  // boundaries
  {register int ind1;
  register int ind2;
  for (ind2 = n2i; ind2 < n2f; ind2++){
      
      // left one
      ind1 = n1r-4;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ( -(*field).u_comp[IND_2(ind1+4,ind2)] 
						   + 4.0*(*field).u_comp[IND_2(ind1+3,ind2)]  
						   - 6.0*(*field).u_comp[IND_2(ind1+2,ind2)] 
						   + 4.0*(*field).u_comp[IND_2(ind1+1,ind2)] 
						   - (*field).u_comp[IND_2(ind1+0,ind2)])*p1/sigma_sphere_rm4;
      ind1 = n1r-3;
      (*dfield).u_comp[IND_2(ind1,ind2)] =   ( 2.0*(*field).u_comp[IND_2(ind1+3,ind2)] 
						    - 9.0*(*field).u_comp[IND_2(ind1+2,ind2)]  
						    + 15.0*(*field).u_comp[IND_2(ind1+1,ind2)] 
						    - 11.0*(*field).u_comp[IND_2(ind1+0,ind2)] 
						    + 3.0*(*field).u_comp[IND_2(ind1-1,ind2)])*p1/sigma_sphere_rm3;
      
      ind1 = n1r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] =   ( -(*field).u_comp[IND_2(ind1+3,ind2)] 
						    + 3.0*(*field).u_comp[IND_2(ind1+2,ind2)]  
						    - 8.0*(*field).u_comp[IND_2(ind1+0,ind2)] 
						    + 9.0*(*field).u_comp[IND_2(ind1-1,ind2)] 
						    - 3.0*(*field).u_comp[IND_2(ind1-2,ind2)])*p1/sigma_sphere_rm2;
      
      ind1 = n1r-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ( - 1.0*(*field).u_comp[IND_2(ind1+3,ind2)] 
						   + 6.0*(*field).u_comp[IND_2(ind1+2,ind2)]  
						   - 14.0*(*field).u_comp[IND_2(ind1+1,ind2)] 
						   + 15.0*(*field).u_comp[IND_2(ind1+0,ind2)] 
						   - 6.0*(*field).u_comp[IND_2(ind1-1,ind2)] 
						   - 1.0*(*field).u_comp[IND_2(ind1-2,ind2)] 
						   + 1.0*(*field).u_comp[IND_2(ind1-3,ind2)])*p1/sigma_sphere_rm1;
      
      // right one
      ind1 = n1s+4;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ( -(*field).u_comp[IND_2(ind1-4,ind2)] 
						   + 4.0*(*field).u_comp[IND_2(ind1-3,ind2)]  
						   - 6.0*(*field).u_comp[IND_2(ind1-2,ind2)] 
						   + 4.0*(*field).u_comp[IND_2(ind1-1,ind2)] 
						   - (*field).u_comp[IND_2(ind1+0,ind2)])*p1/sigma_sphere_sp4;
      ind1 = n1s+3;
      (*dfield).u_comp[IND_2(ind1,ind2)] =   ( 2.0*(*field).u_comp[IND_2(ind1-3,ind2)] 
						    - 9.0*(*field).u_comp[IND_2(ind1-2,ind2)]  
						    + 15.0*(*field).u_comp[IND_2(ind1-1,ind2)] 
						    - 11.0*(*field).u_comp[IND_2(ind1+0,ind2)] 
						    + 3.0*(*field).u_comp[IND_2(ind1+1,ind2)])*p1/sigma_sphere_sp3;
      
      ind1 = n1s+2;
      (*dfield).u_comp[IND_2(ind1,ind2)] =   ( -(*field).u_comp[IND_2(ind1-3,ind2)] 
						    + 3.0*(*field).u_comp[IND_2(ind1-2,ind2)]  
						    - 8.0*(*field).u_comp[IND_2(ind1+0,ind2)] 
						    + 9.0*(*field).u_comp[IND_2(ind1+1,ind2)] 
						    - 3.0*(*field).u_comp[IND_2(ind1+2,ind2)])*p1/sigma_sphere_sp2;
      
      ind1 = n1s+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] =  ( - 1.0*(*field).u_comp[IND_2(ind1-3,ind2)] 
						   + 6.0*(*field).u_comp[IND_2(ind1-2,ind2)]  
						   - 14.0*(*field).u_comp[IND_2(ind1-1,ind2)] 
						   + 15.0*(*field).u_comp[IND_2(ind1+0,ind2)] 
						   - 6.0*(*field).u_comp[IND_2(ind1+1,ind2)] 
						   - 1.0*(*field).u_comp[IND_2(ind1+2,ind2)] 
						   + 1.0*(*field).u_comp[IND_2(ind1+3,ind2)])*p1/sigma_sphere_sp1;
      
    } } 

  //------------------------- add dissipation in the y-direction ----------------------------//
  // interior
  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
    for (ind2 = n2r; ind2 < n2s+1; ind2++){

	(*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	  - ( (*field).u_comp[IND_2(ind1,ind2+4)] 
	      - 8.0*(*field).u_comp[IND_2(ind1,ind2+3)] 
	      + 28.0*(*field).u_comp[IND_2(ind1,ind2+2)]   
	      - 56.0*(*field).u_comp[IND_2(ind1,ind2+1)]   
	      + 70.0*(*field).u_comp[IND_2(ind1,ind2+0)]  
	      - 56.0*(*field).u_comp[IND_2(ind1,ind2-1)]   
	      + 28.0*(*field).u_comp[IND_2(ind1,ind2-2)]
	      - 8.0*(*field).u_comp[IND_2(ind1,ind2-3)]
	      + (*field).u_comp[IND_2(ind1,ind2-4)]
	      )*p2;
	
      
      } } }
  
  // boundaries

  {register int ind1;
  register int ind2;
  for (ind1 = n1i; ind1 < n1f; ind1++){
      
      // left one
      ind2 = n2r-4;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( -(*field).u_comp[IND_2(ind1,ind2+4)] 
	    + 4.0*(*field).u_comp[IND_2(ind1,ind2+3)]  
	    - 6.0*(*field).u_comp[IND_2(ind1,ind2+2)] 
	    + 4.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	    - (*field).u_comp[IND_2(ind1,ind2+0)])*p2/sigma_sphere_rm4;

      ind2 = n2r-3;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( 2.0*(*field).u_comp[IND_2(ind1,ind2+3)] 
	    - 9.0*(*field).u_comp[IND_2(ind1,ind2+2)]  
	    + 15.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	    - 11.0*(*field).u_comp[IND_2(ind1,ind2+0)] 
	    + 3.0*(*field).u_comp[IND_2(ind1,ind2-1)])*p2/sigma_sphere_rm3;
      
      ind2 = n2r-2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( -(*field).u_comp[IND_2(ind1,ind2+3)] 
	    + 3.0*(*field).u_comp[IND_2(ind1,ind2+2)]  
	    - 8.0*(*field).u_comp[IND_2(ind1,ind2+0)] 
	    + 9.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	    - 3.0*(*field).u_comp[IND_2(ind1,ind2-2)])*p2/sigma_sphere_rm2;
      
      ind2 = n2r-1;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( - 1.0*(*field).u_comp[IND_2(ind1,ind2+3)] 
	    + 6.0*(*field).u_comp[IND_2(ind1,ind2+2)]  
	    - 14.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	    + 15.0*(*field).u_comp[IND_2(ind1,ind2+0)] 
	    - 6.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	    - 1.0*(*field).u_comp[IND_2(ind1,ind2-2)] 
	    + 1.0*(*field).u_comp[IND_2(ind1,ind2-3)])*p2/sigma_sphere_rm1;
      
      // right one
      ind2 = n2s+4;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( -(*field).u_comp[IND_2(ind1,ind2-4)] 
	    + 4.0*(*field).u_comp[IND_2(ind1,ind2-3)]  
	    - 6.0*(*field).u_comp[IND_2(ind1,ind2-2)] 
	    + 4.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	    -     (*field).u_comp[IND_2(ind1,ind2+0)])*p2/sigma_sphere_sp4;

      ind2 = n2s+3;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( 2.0*(*field).u_comp[IND_2(ind1,ind2-3)] 
	    - 9.0*(*field).u_comp[IND_2(ind1,ind2-2)]  
	    + 15.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	    - 11.0*(*field).u_comp[IND_2(ind1,ind2+0)] 
	    + 3.0*(*field).u_comp[IND_2(ind1,ind2+1)])*p2/sigma_sphere_sp3;
      
      ind2 = n2s+2;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( -(*field).u_comp[IND_2(ind1,ind2-3)] 
	    + 3.0*(*field).u_comp[IND_2(ind1,ind2-2)]  
	    - 8.0*(*field).u_comp[IND_2(ind1,ind2+0)] 
	    + 9.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	    - 3.0*(*field).u_comp[IND_2(ind1,ind2+2)])*p2/sigma_sphere_sp2;
      
      ind2 = n2s+1;
      (*dfield).u_comp[IND_2(ind1,ind2)] = (*dfield).u_comp[IND_2(ind1,ind2)]
	+ ( - 1.0*(*field).u_comp[IND_2(ind1,ind2-3)] 
	    + 6.0*(*field).u_comp[IND_2(ind1,ind2-2)]  
	    - 14.0*(*field).u_comp[IND_2(ind1,ind2-1)] 
	    + 15.0*(*field).u_comp[IND_2(ind1,ind2+0)] 
	    - 6.0*(*field).u_comp[IND_2(ind1,ind2+1)] 
	    - 1.0*(*field).u_comp[IND_2(ind1,ind2+2)] 
	    + 1.0*(*field).u_comp[IND_2(ind1,ind2+3)])*p2/sigma_sphere_sp1;

    } } 
}
      



/* ------------------------------------------------------------------------------------------ */




/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */


/* ------------------------------------------------------------------------------------------ */

/* Second order dissipation operator like in Carpenter, Nordstrom and Gotttlieb               */

/* NOT FINISHED YET */

void diss_CNG_2d(struct GRID_PAR *grid,
	     struct field *field,
	     struct field *dfield)
{

    int n1i = (*grid).start_grid_1;
    int n1f = (*grid).final_grid_1;
    int n2i = (*grid).start_grid_2;
    int n2f = (*grid).final_grid_2;

    int n_gridpts_1 = n1f-n1i;
    int n_gridpts_2 = n2f-n2i;

    FLOAT xi = (*grid).initial_x;
    FLOAT xf = (*grid).final_x;
    FLOAT yi = (*grid).initial_y;
    FLOAT yf = (*grid).final_y;


    FLOAT d2 =  2.0;

    FLOAT h_1 = (FLOAT)(n1f-n1i-1)/(xf-xi);
    FLOAT h_2 = (FLOAT)(n2f-n2i-1)/(yf-yi);



	/* inner points */
	{register int ind1;
	register int ind2;
	  for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

		(*dfield).u_comp[IND_2(ind1,ind2)] = ((*field).u_comp[IND_2(ind1-1,ind2)]
					    	 - d2*(*field).u_comp[IND_2(ind1,ind2)]
					    	   +  (*field).u_comp[IND_2(ind1+1,ind2)]
						      )*h_1
					           + ((*field).u_comp[IND_2(ind1,ind2-1)]
					    	 - d2*(*field).u_comp[IND_2(ind1,ind2)]
					    	   +  (*field).u_comp[IND_2(ind1,ind2+1)]
						      )*h_2;
		}
		}
	}

 	/* x=0,1 faces */

	{register int ind2;
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

		(*dfield).u_comp[IND_2(n1i,ind2)] =     ((*field).u_comp[IND_2(n1i,ind2)]
						    - d2*(*field).u_comp[IND_2(n1i+1,ind2)]
					    	       + (*field).u_comp[IND_2(n1i+2,ind2)]
							 )*h_1

		                                      + ((*field).u_comp[IND_2(n1i,ind2-1)]
					            - d2*(*field).u_comp[IND_2(n1i,ind2)]
					    	      +  (*field).u_comp[IND_2(n1i,ind2+1)]
							 )*h_2;


		(*dfield).u_comp[IND_2(n1f-1,ind2)] =  ((*field).u_comp[IND_2(n1f-3,ind2)]
						   - d2*(*field).u_comp[IND_2(n1f-2,ind2)]
					    	      + (*field).u_comp[IND_2(n1f-1,ind2)]
							   )*h_1

						     + ((*field).u_comp[IND_2(n1f-1,ind2-1)]
					           - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
					       	     +  (*field).u_comp[IND_2(n1f-1,ind2+1)]
							)*h_2;
		}
	}

	 	/* y=0,1 faces */

	{register int ind1;
	    for (ind1 = n1i+1; ind1 < n1f-1; ind1++){

		(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
					    	- d2*(*field).u_comp[IND_2(ind1,n2i)]
					          +  (*field).u_comp[IND_2(ind1+1,n2i)]
						     )*h_1

						  + ((*field).u_comp[IND_2(ind1,n2i)]
						- d2*(*field).u_comp[IND_2(ind1,n2i+1)]
					    	  +  (*field).u_comp[IND_2(ind1,n2i+2)]
						     )*h_2;



		(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
					    	  - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
					    	    +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
						       )*h_1

						    + ((*field).u_comp[IND_2(ind1,n2f-3)]
						  - d2*(*field).u_comp[IND_2(ind1,n2f-2)]
					    	    +  (*field).u_comp[IND_2(ind1,n2f-1)]
					               )*h_2;
		}
	}




 /* edges */
 
 /* runs on x */

	{register int ind1;
	    for (ind1 = n1i+1; ind1 < n1f-1; ind1++){
(*dfield).u_comp[IND_2(ind1,n2i)] = ((*field).u_comp[IND_2(ind1-1,n2i)]
			       - d2*(*field).u_comp[IND_2(ind1,n2i)]
				 +  (*field).u_comp[IND_2(ind1+1,n2i)]
				    )*h_1

			         + ((*field).u_comp[IND_2(ind1,n2i)]
			       - d2*(*field).u_comp[IND_2(ind1,n2i+1)]
				 +  (*field).u_comp[IND_2(ind1,n2i+2)]
				    )*h_2;


(*dfield).u_comp[IND_2(ind1,n2f-1)] = ((*field).u_comp[IND_2(ind1-1,n2f-1)]
				 - d2*(*field).u_comp[IND_2(ind1,n2f-1)]
				   +  (*field).u_comp[IND_2(ind1+1,n2f-1)]
				      )*h_1

                                   + ((*field).u_comp[IND_2(ind1,n2f-3)]
				 - d2*(*field).u_comp[IND_2(ind1,n2f-2)]
				   +  (*field).u_comp[IND_2(ind1,n2f-1)]
				      )*h_2;
	    }
	}



/* runs on y */

	{register int ind2;
	    for (ind2 = n2i+1; ind2 < n2f-1; ind2++){

(*dfield).u_comp[IND_2(n1i,ind2)] =    (
				       (*field).u_comp[IND_2(n1i,ind2)]
				  - d2*(*field).u_comp[IND_2(n1i+1,ind2)]
				  +    (*field).u_comp[IND_2(n1i+2,ind2)]
				       )*h_1

				  +   (
				       (*field).u_comp[IND_2(n1i,ind2-1)]
				  - d2*(*field).u_comp[IND_2(n1i,ind2)]
				  +    (*field).u_comp[IND_2(n1i,ind2+1)]
				       )*h_2;


(*dfield).u_comp[IND_2(n1f-1,ind2)] =  (
				       (*field).u_comp[IND_2(n1f-3,ind2)]
				  - d2*(*field).u_comp[IND_2(n1f-2,ind2)]
				  +    (*field).u_comp[IND_2(n1f-1,ind2)]
				       )*h_1

				  +   (
				       (*field).u_comp[IND_2(n1f-1,ind2-1)]
				  - d2*(*field).u_comp[IND_2(n1f-1,ind2)]
				  +    (*field).u_comp[IND_2(n1f-1,ind2+1)]
				       )*h_2;
						}
	}



/* corners */

(*dfield).u_comp[IND_2(n1i,n2i)] = (
					(*field).u_comp[IND_2(n1i,n2i)]  
				   - d2*(*field).u_comp[IND_2(n1i+1,n2i)]
				   +    (*field).u_comp[IND_2(n1i+2,n2i)]
					)*h_1

                                     + (
					(*field).u_comp[IND_2(n1i,n2i)]
				   - d2*(*field).u_comp[IND_2(n1i,n2i+1)]
				   +    (*field).u_comp[IND_2(n1i,n2i+2)]
					)*h_2;


(*dfield).u_comp[IND_2(n1i,n2f-1)] = (
					(*field).u_comp[IND_2(n1i,n2f-1)]  
				   - d2*(*field).u_comp[IND_2(n1i+1,n2f-1)]
				   +    (*field).u_comp[IND_2(n1i+2,n2f-1)]
					)*h_1

                                     + (
					(*field).u_comp[IND_2(n1i,n2f-3)]
				   - d2*(*field).u_comp[IND_2(n1i,n2f-2)]
				   +    (*field).u_comp[IND_2(n1i,n2f-1)]
					)*h_2;



(*dfield).u_comp[IND_2(n1f-1,n2i)] = (
					(*field).u_comp[IND_2(n1f-3,n2i)]  
				   - d2*(*field).u_comp[IND_2(n1f-2,n2i)]
				   +    (*field).u_comp[IND_2(n1f-1,n2i)]
					)*h_1

                                     + (
					(*field).u_comp[IND_2(n1f-1,n2i)]
				   - d2*(*field).u_comp[IND_2(n1f-1,n2i+1)]
				   +    (*field).u_comp[IND_2(n1f-1,n2i+2)]
					)*h_2;


(*dfield).u_comp[IND_2(n1f-1,n2f-1)] = (
					(*field).u_comp[IND_2(n1f-3,n2f-1)]  
				   - d2*(*field).u_comp[IND_2(n1f-2,n2f-1)]
				   +    (*field).u_comp[IND_2(n1f-1,n2f-1)]
					)*h_1

                                     + (
					(*field).u_comp[IND_2(n1f-1,n2f-3)]
				   - d2*(*field).u_comp[IND_2(n1f-1,n2f-2)]
				   +    (*field).u_comp[IND_2(n1f-1,n2f-1)]
					)*h_2;


}

#endif /* DISSIPATION */







#endif
